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SPACE-FREQUENCY  CORRELATIONS  IN  MULTISTATIC  ACOUSTIC 
REVERBERATION  DUE  TO  A  WIND-DRIVEN  SEA  SURFACE: 
THEORETICAL  RESULTS  AT  LOW  FREQUENCY 


1.  INTRODUCTION 
1.1  Overview 

For  decades,  the  U.S.  Navy’s  operational  needs  were  well  served  by  the  pursuit  of  acoustic 
antisubmarine  warfare  (ASW)  as  a  largely  passive  enterprise.  In  recent  years,  however,  the  emphasis 
has  shifted  strongly  toward  active  operations,  and  this  change  has  lent  newfound  importance  to 
some  hitherto  secondary  topics.  Prominent  among  these  is  the  irregular  and  dynamic  nature  of 
the  sea  surface.  Previously  of  interest  only  in  connection  with  ambient  noise  and  propagation  loss, 
the  surface  has  taken  on  new  importance  as  a  source  of  reverberation  in  active  systems.  Features 
of  particular  interest  in  the  reverberant  field  are  (a)  Doppler  shifts,  which  can  mimic  or  mask  the 
echoes  from  moving  targets,  and  (b)  two-point  spatial  correlations,  whose  nature  may  ultimately 
help  systems  distinguish  target  returns  from  surface  reverberation. 

This  report  treats  the  second-order  space-frequency  statistics  of  reverberation  from  time- 
dependent  sea  surfaces.  The  focus  is  on  low  frequencies  and  wind  speeds,  where  interface  scat¬ 
tering  dominates  bubble  effects  [1].  Theoretical  techniques  are  used  to  express  the  reverberation’s 
Doppler  characteristics  and  spatial  correlations  in  terms  of  accepted  empirical  parameterizations 
for  the  surface-wave  statistics.  The  overall  goal  is  a  comprehensive  expression  for  the  second- 
order  spatial  and  temporal  correlation  features  that  relate  to  acoustic  systems  operating  in  the 
low-frequency  (LF)  band.  Naturally,  some  judicious  physical  approximations  must  be  made.  The 
source  is  treated  as  a  motionless  point  that  emits  a  narrowband  signal.  The  environment  is  allowed 
to  have  depth  dependence  (including  a  layered  bottom),  but  every  part  of  it  is  range-independent, 
deterministic,  and  static — with  one  exception.  That  exception  is  the  sea  surface,  which  is  composed 
of  waves  traveling  in  all  directions  with  various  (and  ultimately  random)  amplitudes  and  phases.  A 
comprehensive  new  formulation  is  derived  for  the  cross-spectral  density  (CSD)  of  the  reverberation 
in  terms  of  (a)  the  second-order  statistics  of  the  sea  surface  and  (b)  the  Green’s  function  for  the 
same  environment  without  surface  waves.  This  is  expected  to  provide  a  basis  for  future  work  in 
complex  environments;  however,  for  the  remainder  of  this  report,  the  subsurface  medium  is  repre¬ 
sented  as  a  homogeneous  half-space.  This  simplification  allows  analytic  methods  to  reduce  a  crucial 
fivefold  spatial  integration  to  a  single  azimuth  integral.  The  result  is  a  novel  analytic  expression 
from  which  model  computations  are  easily  done.  Sample  computations  are  carried  out  to  display 
the  dependence  of  the  amplitudes  and  phases  of  the  CSD  sidebands  on  the  multistatic  horizontal 
and  vertical  placement  of  the  receiver  pair,  the  wind  speed,  and  the  directionality  of  the  surface 
waves. 

Section  2  deals  on  a  deterministic  level  with  the  two  fundamental  elements  of  the  problem — 
the  sea  surface  and  the  acoustic  field  scattered  from  it.  First,  the  makeup  of  a  deterministic 
time-dependent  surface  is  discussed  and  then  the  small-waveheight  approximation  (SWHA)  is  in¬ 
voked  to  characterize  the  narrowband  signal  that  results  from  scattering  by  such  a  surface.  In 
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Section  3,  the  surface  elevation  is  treated  as  a  random  process  that  is  stationary  both  in  time  and 
in  horizontal  coordinates.  A  SWHA  series  is  obtained  for  the  resulting  reverberation  CSD  through¬ 
out  the  water  and  sediment,  and  the  terms  are  grouped  into  baseband  contributions  (zeroth  order 
and  second  order)  and  Bragg-Doppler  sidebands  (second  order  only).  The  sideband  expression, 
Eq.  (30),  is  recognized  as  an  active-scattering  extension  of  the  van  Cittert-Zernike  theorem  from 
the  classical  theory  of  partially  coherent  fields.  With  the  sea  surface  described  in  terms  of  its  power 
spectrum  and  directional  spectrum,  each  CSD  term  is  reduced  to  an  azimuthal  integration  involving 
a  fundamental  integral  over  surface  wavenumber.  Section  3.4  summarizes  the  results.  In  Section 
4,  the  fundamental  integrals  for  the  baseband  and  sidebands  are  evaluated  by  two-dimensional 
(2-D)  stationary  phase  estimation  for  a  uniform  ocean.  This  uses  a  source/receiver  geometry  that 
yields  analytic  results:  the  source  and  receivers  are  deployed  at  depths  much  less  than  the  ranges 
involved.  Section  5  presents  a  set  of  computer  simulations  for  the  sidebands.  These  are  done  using 
the  Pierson-Moskowitz  power  spectrum  and  Longuet- Higgins  directional  spectrum  to  illustrate  the 
sensitivity  of  the  sidebands  to  receiver  placement,  frequency,  and  sea-surface  parameters.  Section 
6  presents  a  summary  and  conclusions. 

1.2  Background 

Perturbative  treatment  of  sea-surface  scattering  is  certainly  not  a  new  idea.  Indeed,  various 
perturbative  approaches  are  to  be  found  throughout  the  acoustics  and  radar  literatures.  Some  of 
these  deserve  particular  mention  here  in  relation  to  the  present  work. 

In  the  1970s,  Harper  and  Labianca  published  a  series  of  articles  [2-6]  that  applied  analytic 
techniques  to  the  theory  of  bistatic  sea-surface  scattering.  Their  main  result  was  a  formalism  for 
describing  the  Doppler  sidebands  in  the  power  spectral  density  (PSD).  The  present  report  has  cer¬ 
tain  basic  features  in  common  with  that  work.  Both  use  perturbation  theory  along  with  Rayleigh’s 
approximation  for  the  surface  boundary  condition,  and  both  ultimately  treat  the  source  signal  and 
sea  surface  as  stochastic  inputs  that  generate  a  stochastic  output — the  scattered  field — whose  sec¬ 
ond  moment  is  the  basic  goal  of  the  analysis.  In  other  respects,  however,  the  two  approaches  differ. 
Harper  and  Labianca  began  by  modeling  the  input  stochastic  processes  themselves.  They  then 
used  analytic  wave-propagation  techniques  to  determine,  the  resulting  output  process  and,  finally, 
ensemble- averaged  a  quadratic  expression  in  the  output  to  produce  the  desired  second  moment. 
This  required  devoting  considerable  effort  to  devising  an  ensemble  of  faithful  realizations  of  the 
sea  surface  and  to  obtaining  the  consequent  realizations  for  the  perturbative  contributions  to  the 
scattered  field,  including  those  that  ultimately  contribute  nothing  to  the  ensemble  average. 

In  contrast,  this  report  (a)  addresses  multistatic  sea-surface  scattering  and  (b)  approaches  the 
subject  directly  on  the  level  of  second-order  statistics  [7,  Chapter  5].  It  follows  the  basic  agenda  of 
the  classical  theory  of  partially  coherent  wave  fields  [8-11],  exploiting  an  input/output  relation  for 
the  field’s  second  moment  rather  than  for  the  field  itself.  That  relation  arises  as  follows.  The  field’s 
general  space-time  second  moment  r(fi,<i,f2>*2)  obeys  two  wave  equations — in  r\,t\  and  in  r%,t2. 
For  a  field  that  is  statistically  stationary  in  time,  this  moment  reduces  to  the  Mutual  Coherence 
Function  (MCF),  C(ri,r2,  r),  which  depends  on  the  time  difference  r  =  t%  —  <i.  Its  transform  in 
the  frequency  domain  is  the  CSD,  C(fi,  rq,  /),  which  obeys  a  pair  of  Helmholtz  equations — in  rq 
and  in  rq — within  any  source-free  region  r\.  fq  G  V.  Green’s  theorem  then  provides  the  desired 
input/output  relation — an  expression  for  the  CSD  inside  V  in  the  form  of  an  integral  over  the 
boundary  /  q .  rq  £  dV  involving  the  appropriate  Green’s  function  and  boundary  values. 

In  this  form,  that  relation  would  be  prohibitively  difficult  to  use,  since  its  input  involves  in¬ 
tegration  over  the  actual  sea  surface.  Fortunately,  in  the  small-waveheight  regime,  one  can  use  a 
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classic  Rayleigh  technique  to  replace  the  exact  expression  by  a  perturbative  series  approximation, 
each  term  of  which  involves  integration  over  the  mean  sea-surface — a  simple  horizontal  plane.  In 
this  perturbative  representation,  the  nth-order  contribution  to  the  output  is  generated  by  the  nth 
angle  and  frequency  moments  of  the  surface-wave  distribution.  Harper  and  Labianca  followed  ba¬ 
sically  the  same  procedure  in  producing  a  perturbative  series  for  the  field  amplitude.  However,  the 
desired  result — the  nth  perturbative  contribution  to  the  output  moment — is  more  readily  obtained 
when  the  procedure  is  applied  to  the  moments  rather  than  to  the  stochastic  process  itself.  In  that 
context,  it  is  immediately  clear,  for  example,  that  the  first-order  contribution  to  the  CSD  vanishes 
identically.  Working  on  the  level  of  moments  also  seems  to  have  an  advantage  in  terms  of  reliability. 
Their  work  with  the  stochastic  field  amplitude  led  Harper  and  Labianca  to  assert  that  reciprocity 
failed  beyond  the  second  perturbative  order  [6,  p.  1,149].  This  claim  appears  to  be  an  artifact  of 
the  complexity  of  their  approach  and  is  not  borne  out  here. 

In  1993,  Goalwin  published  an  article  [12]  that  modeled  the  multistatic  reverberation  CSD 
for  a  pair  of  receivers.  His  work,  like  the  present  effort,  approached  the  problem  on  the  level 
of  second  moments  and  obtained  a  multiple-integral  expression  for  the  second-order  perturbative 
contribution.  To  facilitate  the  analysis,  however,  he  placed  the  receivers  at  a  common  depth  and 
was,  thus,  unable  to  examine  the  vertical  aspects  of  spatial  correlation.  More  significantly,  he 
dealt  with  the  sea  surface  as  a  “frozen”  (motionless)  interface,  thereby  precluding  any  treatment 
of  surface  Doppler.  In  addition,  like  Harper  and  Labianca,  Goalwin  applied  the  stationary-phase 
method  only  numerically,  obtaining  no  analytic  results  comparable  to  those  produced  here. 

1.3  Conventions  and  Notation 

Vectors  in  3-D  are  indicated  by  arrows,  and  2-D  vectors  and  matrices  are  denoted  by  underlining 
(e.g.,  3-D  vector  r  =  (r,  z)  and  2-D  matrix  ^).  Directions  in  the  horizontal  plane  are  often  specified 
in  terms  of  the  a  unit  vector  n{6)  that  points  along  an  azimuth  9  relative  to  the  dominant  wind 
direction.  The  frequency  /  or  its  angular  counterpart  u  =  2ivf  (or  both)  may  appear  in  any  given 
expression;  the  choice  is  entirely  a  matter  of  convenience.  The  same  is  true  for  the  horizontal 
wavevector  k  =  27 rs. 

The  standard  Fourier  transformation  connects  the  time  and  frequency  domains:  a(f )  = 
f-tt  dt  E(+ft)a(t)  with  the  kernel  E(q )  =  exp(i27rg).  The  same  basic  symbol  is  used  in  both 
domains — a(t),  a(/) — with  only  the  argument  identifying  which  one  is  meant.  The  relation  is  often 

f  f 

denoted  aft )  — 4  a(/),  with  the  arrow  indicating  the  forward  direction.  The  variables  above  the 
arrow  are  sometimes  omitted  when  they  should  be  clear  from  context.  The  spatial  Fourier  trans¬ 
form  in  two  dimensions  follows  the  same  pattern  except  that,  because  of  the  traveling-wave  form 
E(s-r  —  /<),  its  kernel  uses  the  opposite  sign:  a(s)  =  /  d2r_E(-s-r)a(r ).  Again,  the  transformation 
is  denoted  a(r)  =-4  a (5),  with  the  arrow  indicating  the  forward  direction  (this  time,  the  one  with 
the  minus  in  its  kernel). 

When  aft)  is  effectively  constant,  it  will  be  written  as  aft)  to  emphasize  that  it  is  virtually 
a  dc  signal.  To  the  same  approximation,  its  spectrum  is  singular:  a(f)  &  a(/~0)<S(/),  where 
a(/~0)  =  fo-  dfa(f)  denotes  the  integrated  strength  of  the  singularity.  Numerically,  of  course, 
a(/^0)  =  a(jt).  Likewise,  when  afr)  is  effectively  independent  of  r,  it  is  denoted  aff)  and  its 
wavenumber  spectrum  is  a(s)  &  a(s~0)6(s)}  where  a(s~0)  =  aff).  These  are  understood  in 
the  same  spirit  as  expressions  like  “a  <C  /?”  (i.e.,  as  a  physical  notation  rather  than  a  purely 
mathematical  one). 

A  field  a(r,  t)  is  a  function  of  both  spatial  position  and  time.  Double  Fourier  transforma- 
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t  f  V  S 

tion,  a(r,  t)  — 4  a(s,  /),  allows  it  to  be  viewed  as  a  function  of  wavenumber  and  frequency 

a(s.  /).  Although  any  physical  field  can  be  represented  as  a  real- valued  function  in  the  space/time 
domain — a(r,t )  =  a*(r,t),  with  the  corresponding  symmetry  a(s,  f )  =  a*(—s,  —f)  in  the  wavenum¬ 
ber/frequency  domain — it  will  usually  be  more  convenient  to  use  representations  involving  complex¬ 
valued  fields  a(r,t).  This  practice  (the  equivalent  of  representing  cos  q  by  the  real  part  of  etq)  is 
perfectly  harmless  provided  due  care  is  used  with  nonlinear  functions  of  the  field,  particularly  the 
quadratic  form  a*(r1,ti)a(r2,t2)- 

We  will  need  to  deal  with  the  second-order  statistical  moments  of  fields  (i.e.,  ensemble  averages 
of  such  quadratic  forms).  Of  course,  we  may  equally  well  deal  with  moments  of  a(r,  /),  a(s,t),  or 
a(s,f),  instead.  Our  interest,  then,  will  be  focused  on  the  two-point  correlations  in  one  of  these 
forms: 


Ta(r1,ti,r2,t2)  =  (a*(r1,tx)a(r2,t2)) 

r«(£i,  h,r2,  h)  =  (a*(ti,  h)a(r2,  /2)> 

ra(li,  fi,s2 ,  f2)  =  {a*(s i,  fi)a(s2,  f2 ))  . 

Here  too,  the  same  main  symbol  is  used  universally,  and  only  the  arguments  distinguish  the 
space/ wavenumber  and  time/frequency  domains.  Naturally,  since  Fourier  transformations  con¬ 
nect  these  domains,  all  three  forms  are  equivalent.  In  dealing  with  statistically  stationary  fields, 
it  is  often  convenient  to  use  “mean  and  difference”  coordinates  (i.e.,  to  re-express  the  dependence 
on  a  coordinate  pair  in  terms  of  a  dependence  on  their  mean  (ij  —  ^((j  +  (;)  and  difference 

Qj  =  Q  —  £,•).  In  this  representation,  the  three  moments  above  are 

ra(ni2><12,£i2><12)  =  (a*  (rx,t{)a{r_2,t2)) 

ra(l.l2,  /l2,£l2,  /12)  =  /l)«(zi2,  /2)) 

ro(li2,/l2,ll2)/l2)  =  («*(£l,/l)aU2>/2))  • 

Temporal  stationarity  is  so  common  that  special  notation  and  terminology  are  used. 

£<*(2:1,12,  *12)  d-  ^(2:12,  <12,  £12,^12)  =  £0(2^,211, -<12) 

is  called  the  MCF,  and  its  tx— -&2  transform, 

£a(zLi,2l2,/l2)  d-  ra(r12,/i2,£i2,/l2~0)  =  Oz^Ei,  /12), 
is  known  as  the  CSD  or  Mutual  Spectral  Density  Function  (MSDF). 

2.  AMPLITUDES 

We  begin  by  formulating  the  problem  in  the  context  of  a  time-dependent  linear  system  with  one 
input,  e(t),  the  signal  emitted  by  the  source;  one  internal  time  dependence,  h(r,t),  the  sea-surface 
elevation;  and  one  output,  A(r,  t),  the  modulation  of  the  field  scattered  by  the  surface.  Initially 
these  are  all  treated  as  deterministic  quantities.  A  stochastic  treatment  follows  in  Section  3. 

2.1  Surface  Elevation 

Even  if  the  surface  elevation  were  completely  arbitrary,  it  could  still  be  Fourier  analyzed  via 
h(r,t )  -  4^— I  h(s,  f)  (i.e.,  synthesized  from  unrestricted  sinusoids  E(s-r  —  ft))  according  to 

r  n  r+oo 

h(r,t)=  d2s  dfE(s-r-ft)h(s,f),  (1) 

Jr2  J—oo 
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with  only  the  requirement  that  /i(s,  f)  =  /&*(—«,—/)  to  guarantee  that  /i(r,  t)  is  real.  But  of 
course  the  surface  elevation  is  not  that  arbitrary.  Its  components  are  traveling  waves1  that  obey 
the  appropriate  dispersion  relation  for  waves  on  the  sea  surface — a  constraint  in  the  form  of  a 
mapping  |/|  =  F(s),  where  5  =  |s|  (or,  equivalently,  s  —  5(|/|)  in  terms  of  the  inverse  map 
S  =  F^1).  The  functional  form  of  F(s)  depends  on  the  nature  of  the  restoring  forces  acting  on 
surface  displacements  at  the  fundamental  length  scale  of  the  problem.  For  the  present,  we  only 
note  that  F(s )  must  be  an  increasing  differentiable  function  on  s  >  0  and  vanish  at  the  origin 
(F'(s)  >  0,  F(0)  =  0).  This  is  general  enough  to  span  the  spectrum  from  capillary  to  gravity  waves 
in  deep  or  shallow  water  [14,  Eq.  (6.3.12)].  Incorporating  a  surface  dispersion  relation  imposes  the 
restriction  s  =  S(\f\)n(6)  on  the  3-D  integration  region  {s  E  i?2,  |/|  <  oo}  in  Eq.  (1),  shrinking  it 
to  the  2-D  manifold  {s  E  i?2,  |/|  =  F(s)}. 

Alternatively,  one  may  use  a  synthesis  like 

h(r,t)  =  j>  dd  dfa(6,  /)  cos  {<f>(6,f)  +  2w[S(f)n(0)-r-  ft]}  (2) 

in  which  §  d6  is  a  360°  azimuth  integration,  and  the  traveling-wave  component  for  each  {0,  /}  has 
an  amplitude  a  and  a  phase  </>,  both  real.  This  arguably  has  more  down-to-earth  physical  realism 
on  its  side  because  (a)  it  incorporates  the  dispersion  relation  in  a  natural  way,  (b)  it  is  manifestly 
real-valued,  and  (c)  it  avoids  double-counting  (since  /  >  0,  the  component  with  {0,  |/|}  is  included, 
but  its  physical  twin  with  {0  +  7r,  — 1/|)  *s  n°t)- 

Used  consistently,  the  two  approaches  are,  naturally,  equivalent.  The  a  and  <f>  in  Eq.  (2)  are 
specified  only  for  positive  frequencies,  but  their  definitions  can  be  extended  to  the  negative  spectrum 
via  a(0 ,  — /)  =  +a(6  +  tt, /)  and  </>($,  —f)  =  —<f>(Q  +  7 r,/).  In  terms  of  the  complex  amplitude, 
u(0 ,  /)  =  ^a(0,  whose  symmetry  is  u(6 ,  /)  =  u*(6  +  7r,  — /),  Eq.  (2)  is  simply 

/+oo  r 

df  j  d9u(6,f)E(S(f)n(6)-r-  ft)  .  (3) 

Any  physically  realizable  sea  surface  can  be  synthesized  in  this  way  by  integrations  over  azimuth 
and  frequency.  The  symmetry  of  u  guarantees  h(s_,  /)  =  /&*(—£,—/),  so  that  /i(r,  t)  remains  real. 
A  comparison  of  Eqs.  (1)  and  (3)  shows  that 

h(sj)  =  6{s~s(f))u(ej) ,  (4) 

which  illustrates  how  the  dispersion  relation  restricts  the  surface  elevation  in  the  wavenum¬ 
ber/frequency  domain. 

2.2  Scattered  Field 

Initially  we  regard  the  sea  surface,  2  =  /i(r,  t),  as  a  deterministic  boundary  and  present  an 
expression  for  the  acoustic  field  that  scatters  from  it.  The  subsurface  acoustic  environment  is 
prescribed  by  the  sound  speed  c  and  density  p,  which  are  considered  to  be  continuous  functions  of 
depth  except  for  simple  discontinuities  wherever  the  material  properties  change  abruptly  (e.g.,  at 
the  sea  floor  and  at  sediment  interfaces). 


1  We  are  representing  these  surface  components  as  perfect  sinusoids.  Their  shape  is  actually  slightly  distorted  by  the 
orbital  motion  of  the  near-surface  water,  but  the  effect  is  negligible  for  small  waveheights  [13]. 
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2.2.1  Differential  Problem 


This  corresponds  to  a  linear  differential  boundary  value  problem  of  the  familiar  form 

WP(r,t)  =Q(r,t )  •••  z  >  h{r,t ) 

(5a) 

o 

II 

HC 

■  *  ^  =  h(r, ,  t ) 

(5b) 

P(f,t)~*  0  •' 

■  •  |r|  — »  oo 

(5c) 

featuring  a  wave  operator  [15],  W  =  —{p{f)d/dr}-{p~^{r)d/df}-\-c~‘i{r)d‘1/dti,  determined  by  the 
medium’s  sound  speed  and  density.  In  signal  engineering  terms,  the  environment,  in  conducting 
the  signal  from  source  to  receiver,  acts  as  a  linear  filter.  Two  things  drive  the  system:  the  source 
function  Q{r,t )  and  the  sea-surface  elevation  h(r,t).  It  is  the  latter  that  makes  the  problem 
challenging  by  causing  the  ocean  to  act  as  a  time-varying  filter.  Other  conceivable  contributors, 
such  as  time-dependence  or  randomness  in  c  and  p,  are  disallowed  in  the  present  formulation. 

Although  this  problem  is  well  posed  in  a  mathematical  sense  [16],  it  is  appropriate  to  make  at 
least  passing  mention  of  some  of  its  physical  shortcomings.  First,  it  is  implicit  that  c  and  p  must 
actually  be  uniform  near  the  surface  if  there  is  to  be  no  time-dependence  within  the  medium.  This 
means  excluding  all  effects  associated  with  near-surface  hydrodynamics,  such  as  the  advection  of 
air  bubbles  by  the  subsurface  water  [17].  In  fact,  scattering  by  near-surface  bubbles  [18],  whether 
advecting  or  not,  is  omitted  altogether.  Furthermore,  a  more  realistic  surface  boundary  condition 
might  even  be  applied.  An  air/water  impedance  condition,  for  instance,  would  allow  transmission 
into  the  air  and  that  may  actually  be  a  significant  loss  mechanism  in  shallow  water  [19]. 

2.2.2  Narrowband  Emission  from  a  Point  Source 

We  suppose  that  the  source  is  fixed  in  position  and  that  its  time-dependence  is  uniform 
throughout.  This  time-dependence  we  take  to  be  a  modulation  applied  to  a  carrier  frequency 
/0.  Thus  the  source  term  is  Q(r,  t )  =  q(r)e(t)  exp(— The  resulting  acoustic  field  has  the  form 
P(r,t )  =  A(r,  t)  exp(-iu)Qt) ,  in  which  A(r,  t),  like  eft),  is  a  slow  complex  modulation.2  We  further 
idealize  the  source  as  a  spatial  point  (i.e.,  take  q(r)  =  8{r  -  f0)).  To  zeroth  order  (i.e.,  neglecting 
dA/dfjJot)  and  <92 A/ diw^t)2  relative  to  A),  we  obtain 

WA(r,t)  =  e(t)6(r  — ro)  ■■■  z  >  h(r,t)  (6a) 

A{f, t)  =  0  •••  z  =  h(r_,t)  (6b) 

A(r,  t)  —+  0  •  •  •  |f|  — *  oo  ,  (6c) 

in  which  >V  =  ~({p{r)d  /  dr)-{p~l  (r)d  /  dr] +kl{f))  is  a  general  Helmholtz  operator  with  a  spatially 
varying  wavenumber  ko(r)  =  wo/c(r).  The  time  variation  in  A  is  caused  by  e  and  h. 

When  h  is  frozen,  the  function  A(r,t)/e{t)  loses  all  dependence  on  time  and  is,  in  fact,  simply 
the  space-frequency  Green’s  function  for  that  fixed  boundary  at  the  frequency  fo  (i.e.,  the  frequency 
response  function  of  the  now  time-invariant  environment  at  the  carrier  frequency).  However,  we 
want  to  deal  with  situations  in  which  h  is  time  dependent.  Our  first  step  in  doing  that  is  to  recast 
Eq.  (6)  in  integral  form. 

2.2.3  Integral  Formulation 

For  any  volume  V,  provided  its  boundary  8V  remains  fixed  and  the  density  is  uniform  along 
the  normal  direction  at  the  boundary,  WU(r)  =  6 (r-  fo)  is,  by  Green’s  theorem,  equivalent  to 


2  A  also  depends  parametrically  on  fo  and  fo,  of  course,  but  this  will  be  left  implicit. 
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U(F)  =  G(r, f0)  -  [  d2rb Gb(r, n)U(rb)  ■■■  r,r0  e  V  , 

where  G(ra,r^)  is  the  Green’s  function3  for  pressure-release  conditions  on  dV> 
notation  used  in  the  integrand  means 

Gb(r,rb)  d=  n{rb)-  Jim  . . .  f  e  y  ,  fb  <E  dV  ,  (8) 

U-*Vf,  V  UU  ) 

with  the  limit  taken  from  within  V.  It  is  simply  a  shorthand  for  the  outward  normal  derivative 
with  respect  to  the  subscripted  argument  (i.e.,  the  “boundary  Green’s  function”  [20]).  Reciprocity 
is  embodied  in  the  expression  G(ra,  r&)p(r&)  =  G(r&,  ra)p(ra). 

Equation  (7)  can  be  applied  to  U  =  A/e  by  identifying  V  with  the  half-space  z  >  0  and  dV 
with  the  z  =  0  plane,  S.  The  result  is  that  throughout  the  water  and  sediment, 

A(r,t)  =  e(t)G(r,r0)  -  [  d2raGa(r,ra)A(ra,t )  .  (9) 

J  s 

G  is  the  “sea  state  zero”  Green’s  function  for  pressure-release  conditions  on  5,  so  that 


(7) 


and  the  subscript 


Ga(r,ra)  =  - 


dG(f,fa)\ 

dza  )Za= 0 


is  just  its  upward  derivative  at  the  surface  point  ra  =  (z^50).  Equation  (9)  relates  A  to  its  values 
on  the  fixed  plane  S',  but  it  is  not  a  solution  because  Eq.  (6)  needs  boundary  values  specified  on 
the  actual  moving  sea  surface. 


2.2. 4  Small- Waveheight  Approximation 


To  get  such  a  solution,  one  must  infer  the  values  of  A(r,  t)  at  z  =  0  from  the  vanishing 
boundary  values  imposed  at  z  =  h(r,  t).  This  can  be  done  by  the  SWHA  [21].  Since  this  is  such  a 
well-known  procedure,  we  omit  the  details  and  simply  note  that  the  modulation  emerges  as  a  series 
A  =  +  A ^  +  A H whose  initial  terms  are 


,4(0)(f,£)  =  e(t)G(r,r0)  (10a) 

A{l\r,t)  =  -e{t)  (  rf2raGa(f,fa)/i(r^,<)Ga(fa,f0)  (10b) 

Js 

A(2\r,t)  =  e(t )  ff  d2r^d2riGa(r,ra)h(za,t)Gab(ra,rb)h(rJ),t)Gb(fb,fo )  .  (10c) 

Again,  the  subscript  notation  of  Eq.  (8)  is  used  for  the  surface  normal  derivatives.  All  of  these 
A terms  have  an  implicit  parametric  dependence  on  tq  and  /o,  and  they  all  derive  their  time 
dependence  from  both  e  and  h  (except  which  is  independent  of  h ).  Furthermore,  each  A ^ 
inherits  the  Green’s  function’s  reciprocity  property.  It  has  been  suggested  [6]  that  this  order-by- 
order  reciprocity  might  not  extend  beyond  second  order;  however  reciprocity  has  been  explicitly 
demonstrated  for  all  orders  in  the  cw  (continuous  wave)  case  with  a  frozen  surface  [22],  and  it 
seems  clear  enough  from  the  above  that  it  does,  in  fact,  persist  in  general.  This  is  the  SWHA 
for  the  complex  modulation  that  results  when  a  given  narrowband  signal  scatters  from  a  given 
time-dependent  wave  surface.  It  appears  here  not  as  a  novel  result  (indeed,  Eq.  (10)  is  a  minor 
variation  on  Eqs.  (30)  and  (31)  of  Ref.  6)  but  rather  as  a  natural  starting  point  for  the  following 
section’s  systematic  treatment  of  the  source  signal  and  the  sea  surface  as  stochastic  processes. 

3 The  space/frequency  Green’s  function’s  dependence  on  /0  is  left  implicit. 
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3.  MOMENTS 

For  most  real-world  applications,  one  must  abandon  the  depiction  of  e(t),  h(r_,  t)  and  A(r,t ) 
as  deterministic  quantities  and  instead  treat  them  as  random  processes.  We  do  that  here,  taking 
all  of  them  to  be  statistically  stationary  so  that  their  mean  values  (e),  ( h ),  (A)  can  be  assumed  to 
vanish  and  attention  can  be  focused  on  their  second  moments. 

3.1  Source  Modulation 

The  narrowband  modulation  e(t)  can  be  characterized  by  its  PSD,  Ce(f).  This  can  be  written 
in  the  form  Ce(f )  =  PeRe(f)  with  the  normalization  df  Re{f )  =  1  so  that  the  emitted  power 

is  Pe  =  dfCe(f  ).  The  Gaussian  case  with  bandwidth  A/,  for  example,  would  have  Re(.f)  = 
exp[-^(//A/)2]/(\/27r A/).  But  whatever  the  statistics,  when  the  bandwidth  vanishes,  Re(f) 
reduces  to  its  dc  limit,  6(f).  We  will  often  use  such  an  ideal  reference  source — a  cw  source  with 
unit  power  output  (i.e.,  Re(f)  =  6(f)  and  Pe  =  1). 

3.2  Surface  Elevation 

Because  the  sea  surface  is,  by  assumption,  temporally  and  spatially  stationary,  its  component 
waves  turn  out  to  form  what  is  sometimes  called  a  “free  wave  field”  (i.e.,  a  collection  of  2-D 
plane  waves  that  have  no  frequency  or  angle  correlation  with  each  other  [23]).  Such  a  surface  is 
customarily  described  in  terms  of  its  frequency  spectrum  and  its  directional  spectrum. 

3.2.1  Free-Wave  Sea  Surface 

In  a  stochastic  description,  the  complex  amplitude  u(6,  f)  becomes  a  random  process.  Since 
(u)  vanishes,  the  relevant  statistic  is  the  second  moment 

Tu(o1,h,e2,f2)  =  (u*(6i,fi)u(e2j2)) .  (li) 

From  the  traveling- wave  surface  representation,  Eq.  (3),  we  have 

/P  /*+00  Z’  +  OO 

dO i  (f>  dO 2  /  df\  /  d/2  {u*(0i,fi)u(O2,f2)) 

xE  (-S(fi)n(9i)-r1  +  /i<i  +  S(f2)n(02)-r2  ~  hh) 

=  J d2sid2s2  dfi  df2  E(---)  ~^(/i))  ^(*2  Ipl)r U(0uflt  02,  /2)  ,  (12) 

with  Sj  =  Sjn(0j)  and  (•  •  •)  =  (I12T12  +  *l2‘£l2)  “  (/l2*12  +  /l2?12)  in  the  second  form.  We  are 
supposing  that  h(r,t)  is  statistically  stationary  in  t,  as  is  commonly  held  to  a  good  approximation. 
Consequently,  r^(zii2)  tl2,f.i2>  ^12)  needs  to  be  independent  of  i\2,  and,  thus,  the  E(—fi2i\2)  factor 
in  the  integrand  cannot  actually  introduce  any  t\2  dependence,  no  matter  what  value  /12  may  take. 
In  other  words,  there  must  be  a  6(f\2)  factor  lurking  in  the  moments.  We  are  also  supposing  that 
h(r_,t)  is  stationary  in  r.  By  similar  reasoning,  the  moments  must  also  contain  the  factor  d(s12)  = 
6(0 1  —02)6(si  —  s2)/si.  But  since  we  already  know  that  a  factor  6(si  —  S(fi))6(s2  —  S(f2))6(fi  —  f2) 
is  present,  the  only  news  provided  by  this  spatial  stationarity  is  that  the  moments  must  contain 
an  angular  singularity  6(0i  -  02)  as  well.  In  other  words,  the  traveling-wave  components  of  the 
surface- wave  field  are  uncorrelated  in  both  frequency  and  azimuth  (i.e.,  they  are  a  free  wave  field). 
For  such  a  field,  it  is  not  difficult  to  show,  via  the  polar  form  s12  =  s\2n(dx2),  that 

Ch(h2,  fn)  =  ^(g12  _  lg^12^C'M(0i2,012^O,  /12)  . 

S 12 

The  form  and  content  of  Eq.  (13)  are  a  direct  inheritance  from  Eq.  (4). 


(13) 
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3.2.2  Surface  Frequency  Spectrum  and  Directionality 

It  is  “customary”  [24]  to  represent  the  above  angle-frequency  moment  in  the  factored  form4 


C«(^12)  $12~0,  /12)  =  $(^12) 


H(0 12) 

H  (012  +  7r) 


•  •  •  2  >  0 

•  •  •  Q 12  <  0  ’ 


(14) 


in  terms  of  the  frequency  spectrum  and  directional  spectrum  H ,  both  of  which  are  real  and 
positive.  We  use  this  representation  throughout  the  theoretical  development  without  specifying  any 
functional  forms  for  $  or  H. 

3.3  Scattered  Field 


We  turn  now  to  the  statistics  of  the  reverberant  field.  For  simplicity,  we  hereafter  drop  the 
subscript  “A”  and  write  T(-  •  •),  C( •  •  •)  rather  than  Ta(m  •  •),  Ca{'  *  •)■  The  modulation  A(r,t)}  like 
the  surface  /i(r,  t)  that  scatters  it,  is  a  stochastic  field — a  random  function  of  space  and  time,  albeit 
a  complex- valued  one.  Our  objective  is  to  find  the  second  moments  of  A{f,  t)  for  various  r,  t.  We 
write  the  general  two-point  space/time  moment  as 

<^*(n,<i)j4(f2,<2)>  =  r(ti, <12, £2**12;  £0,20,21,32),  (15) 

consigning  the  source  position  (r^,  zq)  to  the  list  of  parameters  along  with  the  receiver  depths.  The 
brackets  (  )  denote  averaging  over  ensembles  of  both  source  signals  and  sea  surfaces. 

The  SWHA  series  for  A  immediately  yields  a  similar  series  for  the  second  moment.  With  the 
arguments  omitted  for  clarity,  this  is 

r  =  r*°)  +  rw  +  r<2>  +  •  •  •  .  (i6a) 

The  Nth  term  is  =  En>d> o^(N,D\  where  T^)  =  TJ’m  +  with  N  =  j  +  m  and 

D  =  \j  —  m|,  in  terms  of  the  elementary  (j  +  m)th-order  moments, 

Tjm  =  .  (16b) 

The  terms  through  second  order  are  =  T^0,0),  and  =  I^2,0)  +  I^2,2).  Thus, 

r(°)  -  r00  (17a) 

r(1)  =  r01  +  r10  (17b) 

r(2)  _  p11  +  (r02  +  r20) .  (17c) 

Equation  (17)  clearly  holds  in  the  space/frequency  and  wave  numb  er/frequency  domains  as  well. 
The  space/frequency  version  is  obtained  through  the  double  transformation 

r(Ei,t12,r2,fi2;  •••)  tl2-l12  r(d,/i2,r2,/12;  •••) ,  (18) 

which,  in  the  time-stationary  case,  reduces  to 

C{Ti,r2,ti2;  •••)  C(ei,E2,/i2;  •  ■  •)  .  (19) 

The  SWHA  contributions  to  this  CSD  are  investigated  below  for  orders  N  —  0,1,  2.  As  should  be 
expected  from  the  symmetry  of  the  problem,  these  depend  on  only  the  relative  receiver  locations, 


4The  positive-frequency  part  is  Eq.  (3.22)  in  Ref.  24,  where  the  left-hand  side  would  be  written  F(u>  12,^12). 
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22jqi ,  t^q2  rather  than  on  r0,  r1?  and  r2,  separately.  Anticipating  that,  we  hereafter  write  this  CSD  as 
C(r oi,Io2j  /12;  *  *  *)■  In  a  similar  spirit,  we  often  write  the  Green’s  function  G(rJ,  fj)  as  G(r^;  2;,  zy) 
to  emphasize  its  dependence  on  the  relative  horizontal  separation  r_ji-  The  main  focus  is  on  the 
second-order  CSD,  which  consists  of  two  kinds  of  terms:  and  It  turns  out 

that  the  first  type  merely  supplies  a  small  baseband  correction,  while  the  second  provides  something 
quite  different— Doppler  sidebands. 

3.3.1  Order  0 

From  Eq.  (10a),  we  have 

£0,20,21, z2)  =  <^°)*(fi,ti)^°)(f2,<2)>  =  Te(tut2)G*(rur0)G(r\r0)  .  (20) 

The  stationarity  of  e(t )  converts  this  to 

C^(roi,£o2,/;  z0,zi,z2)  =  Ce(f)G*(rm;  21,20)G(r02 ;  z2,z0)  ,  (21) 

which  is  plainly  a  baseband  term  (one  whose  frequency-dependence  reduces  to  6(f)  if  the  source 
bandwidth  vanishes)  consisting  of  the  source  modulation  PSD  modified  by  a  pair  of  cw  propagation 
factors  at  the  carrier  frequency.  As  anticipated,  its  horizontal  dependence  is  on  Zloi  5  Ho2- 

3.3.2  Order  1 

Because  the  source  signal  and  sea  surface  are  statistically  independent,  their  joint  moments 
factor  ( h  e)  =  ( h){e ).  Since  ( h )  =  0,  vanishes  identically. 

3.3.3  Order  2,  Baseband  Terms 

The  two  terms  comprising  iS2,2)  are 

r02(r1,<i,r2,t2;£o>  z0,  Zi,  22)  =  (A^'>*(r1,t1)A^\r2,t2))  (22a) 

r20(r1,ti,r2,t2]r0,  z0,  Zi,  z2)  =  (>l(2)*(fi,  <i)j4(0)(f2,  t2))  •  (22b) 

We  obtain  the  first  one  in  detail  and  then  get  the  second  by  inspection.  From  Eq.  (10), 
r02(£i,ti,£2,t2;  £0,20,21,22)  =  (e*(ti)e(t2)}G*(r1,ro) 

xff  ^2£a^2£iGa(r2,fa)Ga6(ra,r6)G6(f6,f0)  {h{taM)KtbM))  >  (23) 

where  the  two  ensemble  averages  refer  separately  to  the  source  and  the  sea  surface.  We  rewrite  all 
three  moments  in  mean-and-difference  coordinates  and  then  invoke  stationarity  in  time,  space,  or 
both  (indicated  by  ^)  to  simplify  them. 

•  For  the  source  modulation, 

(e*(<i)e(<2)>  =  r«(ti,t2)  =  r.(ti2,«i2)  -  re(ti2,/12)  =  ce(t12)  ce(f ) . 

•  For  the  surface  elevation, 

^a)^(2l7>5  ^b)}  ~  ^h(lia >  "Lbi  ^b)  l'/i(ZLa&?  ^ab->  I-abi  ^ab)  ^h(?Labi 

/oo  _  _  _ 

dfabj>  d0abE  (S(fab)n(eab)-r^b  -  fabtab)  Cu(9ab,  6ab~ 0,  fab )  . 
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•  For  the  reverberation  field, 

r02(r1,ti,r2,t2;  •••)  =  r02(r1,<i2,r2,<i2;  •••)  r02^, <12,12 .A2;  •••) 

=  c,02(r01,r( 02,  <12;  •■•)  tlJ-J  C02(r 01,  £02,/;  •  •  •)  • 

These  three  expressions  convert  Eq.  (23)  into 

C02(rm,rj02,f;zo,zi,z2)  =  Ce(f)G*(fi,f0) 

x  JJ  d2r^d2r^Ga{r2,ra)Gab{ra,rb)Gb{rb,rQ) 

x  I-ocdf'f  d6E  cu{0,eab~0,f).  (24) 

As  anticipated,  this  is  another  baseband  term  like  Eq.  (21),  though  smaller  and  with  a  more 
complicated  spatial  dependence.  With  a  little  rearrangement,  it  can  be  written  as 

C02(r0x,Lo2,f;zo,zi,z2)  =Ce(f)G*(n,rf0)J^df'f  d0Cu(0,6ab~0,  /') 

X  [£(s,  I02;  *0, 22)]1=5(/<)„(0)  ,  (25) 

in  terms  of  the  mean-sea-surface  fundamental  integral, 

L(s,rJ02’,zo,Z2)  =  J  d2rj,Gb(rj0b-zb,zo)J(s,rl2-,Z2)  ,  (26) 

in  which  J  is  the  convolution 

J  fa  r;  z)  =  Gab(s;  za,zb)  %  {F(-s-r)Ga(s;  z,  za)j  .  (27) 

C20(r01,r02,f;z0,z1,Z2)  comes  from  Eq.  (22b)  by  the  same  process.  Details  are  omitted.  Instead, 
both  terms  are  summarized  in  Section  3.4. 

3.3.4  Order  2 ;  Sideband  Terms 

The  only  term  contributing  to  r(2}0)  is  T11.  Thus, 

r(2’0)(r1,ti,ri2,<2;  rJO,zo,zi,z2)  =  T e(h,t2) J j d2r^d2j^Gl(r^1;zi,za)G*a(rj0a;za,z0) 

x  Gb(ri2;z2,zb)Gb{ri)b]zb,zo)rh(r^,ti,ri,t2)  .  (28) 

In  mean-and-difference  time  coordinates  this  can  be  written  as 

r( 2y0\rl,ti2,r2yii2]  ■■•)  =  re(ti2, ii2)CTh(l«, *12, Z^>, *12)  , 

in  which  C  represents  the  integral  operator  and  the  ellipsis  stands  for  the  parameters  r0,  zo>  zh  z2- 
The  stationarity  of  the  source  implies  Te(ti2:ii2)  =  C*e (^12) ?  the  stationarity  of  the  sea  surface  im¬ 
plies  r /j  (zia,  ^12,  Zl^,  ^12)  =  C/i(l«,l6>^12))  and  together  they  imply  that  A^l\r,  t)  is  stationary  in  time 
and  that  r(2’°)(ri,  ti2,r2,?i2;  •  •  •)  =  i,Io2>^12;  *  *  ■)•  I n  the  frequency  domain,  therefore, 

C('2'0)(Lo1,rj02,f;zo,zi,z2)  =  Ce(f )  *  cj^fam r<)2>/;  *D,*1,  *2)  (29) 

in  which  the  source  PSD  is  convolved  with  the  reference  CSD  (Section  3.1), 

def0)  fan > ro2»  /;  *o,  *1 ,  *2)  = 

//^  d  Tjy  Ga{r_a^^  Z\,  Zq)  ^Ga(?loa;  ^0)^/1  (2^65  /)^&(llo&?  ^0)  ^  ^6(z^2i  ^2?  2&)  .  (30) 
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Equation  (29)  is  not  a  baseband  term:  its  frequency-dependence  does  not  reduce  to  6(f)  in  the 
limit  of  a  cw  source. 

For  a  passive  ASW  problem  with  the  sea  surface  acting  as  the  primary  acoustic  source,  the 
curly  brackets  in  the  integrand  above  would  contain  nothing  but  the  boundary  source  CSD  [11, 
Eq.  (3.19)].  Thus  Eq.  (30)  would  reduce  to  a  form  of  the  van  Cittert-Zernike  theorem  from  partial 
coherence  theory  [9,  Section  10.4.2],  which  gives  the  CSD  within  a  volume  in  terms  of  its  values 
on  the  bounding  surface.  This  is  often  a  convenient  starting  point  for  passive  problems  (e.g.,  the 
Kuperman  and  Ingenito  treatment  of  the  ambient  noise  field  generated  by  the  sea  surface  itself  [25, 
Eq.  (7)]).  In  the  active  problem,  the  sea  surface  is  a  secondary  source — a  scatterer.  Rather  than 
that  simple  boundary  value,  the  curly  brackets  contain  (a)  two  Green’s  functions  that  propagate 
the  signal  from  the  primary  source  to  a  pair  of  surface  points  and  (b)  a  wave-height  CSD  that 
embodies  the  sea-surface  correlation.  In  that  sense,  Eq.  (30)  is  a  perturbative  generalization  of  the 
van  Cittert-Zernike  theorem  to  active  scattering  from  a  time-dependent  surface. 

Further  reduction  of  Eq.  (30)  is  possible  when  Ch  in  the  integrand  is  examined  in  the 

wavenumber  domain.  Generally,  C^T-abiLabi  f)  — ~*b  Ch(!abi  $-abi  /)>  however,  when 

h  is  spatially  stationary  (i.e.,  when  Ch(labi^ab- /)  ~  Ch  (!<,&,  f  )6(sclh))  this  reduces  to 

Ch(r^tbJ)  =  Id%bE(sub-rab)Ch(sab,s<lb^0,f).  Then  Eq.  (30)  becomes 

(£01)1102!  /;  zo,  zi,z2)  =  J  d2s^h  M2Ch(lab,  S«6~0,  /) 

in  terms  of  the  surface  integral  M £  =  f  d2r^E(s^b-r_c)Gc(zce:iZ£,zc)Gc(r^c]zc,ZQ).  It  is  a  simple 
matter  to  express  the  latter  as  Mi  =  E(s_ab-r_Q)Nn  in  which  Ni  is  a  spatial  convolution  which 
becomes  a  simple  product  in  the  wavenumber  domain: 


Ne  =  E(scb-rJ0i)Gc(rOe',  zc,  zo )  **  Gc(r 0£;  zt,  zc )  0L^?e  Gc(sm  -  s^b\  zc,  z0)Gc(soe]  ze,  zc )  . 

Thus,  for  any  real  /,  whether  positive  or  negative,  Eq.  (30)  becomes 

C^izQx,^  f\zo,zi,z2)  =  J  d2sabCh(sab,s<lb^O,f)I*(sab,rjai;zo,zi)I(sab,r02-.lzo,Z2 ) 
in  which 

d(s,  r;  zq,  z)  =  J  d2aE(a-r)Ga(a- s]za,zo)Ga(a;z,za)  .  (31) 

Using  Eq.  (13),  the  Ch  factor  can  be  expressed  as  an  angle/frequency  moment  and,  when  that 
moment  has  the  factored  form  of  Eq.  (14),  the  result  is 


C<£f\r01,rj02,±\f\-,zo,zi,Z2)  =  $(|w|)^  d0H(O)  [I*  (s:rQX]  zq,  z\)I(s,rQ2]  zq,  Z2)]±=±S(\f\)n(e)  >  (32) 

an  expression  with  explicit  upshifted  and  downshifted  Doppler  sidebands. 

Equations  (31)  and  (32)  should  prove  computationally  useful  whenever  the  Green’s  function 
is  simple  enough  for  the  I  terms  to  be  evaluated  analytically  since  that  would  leave  only  a  single 
angular  integration  to  be  done  numerically.  This  plan  is  pursued  in  Section  4,  with  the  ocean 
modeled  as  a  uniform  half-space.  But  first,  the  results  are  gathered  and  summarized  in  Section  3.4, 
using  the  familiar  “/c-w”  notation. 


Space- Frequency  Correlations  in  Multistatic  Surface  Reverberation 


13 


3.4  Summary  of  the  Reverberation  CSD 

We  have  found  that  the  reverberation  CSD  assumes  the  form 

f 

C(rm,rJ02,f;zo,zuZ2)  =  Ce(f)  *  (33) 

because  the  sea  surface  and  the  source  are  temporally  stationary.  We  focus  on  Cre f  since  the  PSD  for 
any  narrowband  source  can  always  be  included  later  by  convolution.  From  the  angular  singularity 
in  Ch  (a  consequence  of  the  spatial  stationarity  of  h),  the  angle-frequency  factorization  ansatz  for 
Cu  (Eq.  (14)),  and  the  horizontal  invariance  of  the  subsurface  medium,  we  have  seen  that  Cref  has 
the  structure  outlined  below  through  second  order  in  the  SWHA. 

3.4-1  Baseband 

There  is  a  zeroth-order  baseband  component, 

c£f  (£oi> l02>/; *0, zi,z2)  =  6(f)G*(rJ01;zi,z0)G(r02]Z2,zo)  .  (34) 

There  is  also  a  pair  of  second-order  baseband  contributions: 

C%(L01,L02j;z0,z1,z2)  =  6(f)G*(fur0)  j°°Jf  §  dOCu(6, 0o6~O,  f) 

x  [L(K_,  Loz'i  z0,  Z2)]K_-K(uj')n{e)  (35a) 

C*}(£01,£02,  f;zo,zhZ2)  =  6(f)G(f2,r0)  J™  df'  f  d6Cu(6, 0ab~O,  /') 

x  [L*(K_,  ZLoi ;  z0,  zl)]K_=K(to’)n(d+ir)  >  (35b) 

in  which 

L(K,  Lot;  zo,  z£)  =  f  d2rj,  Gbirjoi,-,  zb,  z0)J (K,  zt)  (36) 

and 

J(K,  r;  z)  =  — ^  /  d2k  e+i(-~—)'LGab(ki  za,  zb)Ga{K  -  k;  z,  za)  .  (37) 

(27T)Z 


3.4-2  Sidebands 

Nonbaseband  contributions  emerge  at  second  order  as  a  pair  of  Doppler  sidebands, 
C^f\r01,rj02,±\f\;z0,z1,z2)=  $(|w|)  ^d0tf(0) 

X  [e -12  — (JC>  Zloi!  zq,  Z\)I (K_,  Tjq2]  Zo>  z%)\ K=±Kn(0) 

•••  \u\=Sl(K),  K  >  0  (38) 

in  which 

I(K,r;,z0,z )  =  j  d2kei^Ga(k-K/2-za,z0)Ga(k  +  K/2-,z,za)  .  (39) 

Since  the  wavevector  in  Eq.  (38)  is  restricted  to  K_  =  +Kn(8 )  and  K_  =  — ATn(0),  each  surface 
wavenumber  K  >  0  contributes  explicitly  to  both  the  upper  (+|u;|)  and  lower  (—  |u;|)  sidebands. 
For  each  K ,  Eq.  (38)  is  evaluated  by  integrating  over  azimuth  and  then  assigning  the  result  to  the 
proper  Doppler-shifted  frequency,  ±[/|. 
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In  Eqs.  (34-39)  we  have  concise  expressions,  valid  through  second  order  in  the  SWHA,  for  the 
baseband  and  sideband  parts  of  the  reverberation  CSD.  The  zeroth-order  baseband  term  does  not 
relate  to  the  sea  state  at  all.  The  second-order  baseband  and  sideband  terms  do  involve  the  sea 
surface,  but  in  different  ways.  The  baseband  term  only  depends  on  the  average  surface  roughness, 
which  is  constant  in  time,  whereas  the  sideband  term  involves  the  full  surface  dispersion  relation. 
(In  photographic  terms,  a  snapshot  could  determine  the  baseband,  but  the  sideband  would  require 
a  movie.)  The  remaining  task  is  to  determine  the  K_  dependence  of  the  fundamental  integrals  L  for 
the  baseband  and  I  for  the  sidebands.  The  following  section  does  this  for  an  ocean  modeled  as  a 
uniform  half-space. 

4.  SHALLOW  DEPLOYMENT  IN  A  UNIFORM  OCEAN 


The  CSD  expressions  obtained  above  will  now  be  evaluated  assuming  a  uniform  ocean 
and  a  geometry  in  which  all  the  depths  are  much  smaller  than  the  source-receiver  ranges 
(zq,  zi,Z2<^roi,ro2) — a  case  that  retains  considerable  operational  relevance  while  remaining  analyti¬ 
cally  tractable.  The  Ith.  receiver’s  horizontal  displacement  from  the  source  is  written  r0p  =  ropn(0p), 
and  its  azimuth  Op  serves  as  a  convenient  reference  for  surface  waves,  dp  =  0  -  Op  (e.g.,  in 
n(0)-n(0p)  =  cos'd  p). 

4.1  Green’s  Function 


For  an  ocean  modelled  as  a  uniform  half-space  beneath  a  pressure-release  plane,  the  Green’s 
function  is  a  sum  of  source  and  image-source  terms: 


gifco  \/r2+(z'-z)2  eik0  i/r2+(z'+z)2 

^  47 r-^/r2  +  ( ~z '  —  ~zp  47ryV2  +  (F  +  z )2 


(40) 


In  the  horizontal  wavenumber  domain,  this  becomes 

eiK(k0,k)\z'-z\  eiK.(ko,k)\z'+z\ 

G(k,z,z)  2*/e(&o,&)  2iK,(ko,k)  ’ 


(41) 


whose  terms  correspond  to  propagation  with  3-D  wavevectors  k  =  k  +  nez  and  k  =  k  —  nez.  These 
wavevectors’  magnitudes  and  those  of  their  horizontal  and  vertical  components  are,  respectively, 
k0  =  |£|  =  wo/co,  k  —  |&|,  and  K,(ko,k)  =  Jk%  —  k2.  The  surface  normal  derivatives  are  simply 
Ga(k]zb,za)  =  e,K(fc°’fc)^,  Gb(k;zb,za)  =  c**(*o »*)*•,  and  Gab(k]zb,za )  =iK,(k0,k). 

4.2  Baseband,  Order  0 

For  shallow  deployment  (z£,  zo^.rop),  the  Green’s  function  becomes 


G(rjop-,  ze,  z0 )  « 


piko  ro£ 

- : - sj[n 

2mr0i 


/  =  1,2  , 


and  this  reduces  the  zeroth-order  baseband 


component  Eq.(34)  to 


(42) 


cS?(£01»£02>  f]Z0,Zl,Z2) 


eiko(ro2~r01)  /  z  \  / 

6(f)-T2 - sin  \k°z° — )  sin  l1 

47r^r0ir02  V  roiJ  \ 


k0zo 


z2 

ro2 


(43) 


As  should  be  expected,  this  depends  on  the  magnitudes  of  roi,  £o2;  but  not  on  their  directions. 
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4.3  Baseband,  Order  2 

4-3.1  Fundamental  Integral 

In  this  environment,  Eq.  (37)  reduces  to 

J(K,r;z )  =  jf—K  [  d2kK(k0,k)ei^-^+<k°^-^ 

(27 r)zJ 

so  that,  through  Eq.  (36)  with  the  integration  variable  shifted  from  to  r^,  L  becomes 
L(K,Loe-,z0,ze)  =  j  d2kK(k0,k)elK('k°’^~Ei)zeJ d2r^eGb(r^e  -  ^  zb,  z0)e~^K~k'>'L^  .  (44) 
The  spatial  integral  in  this  is  e^-~— * '£oeetK(k°’\--—\ D*°,  so  the  fundamental  integral  is  just 

L(K, Loe;  z0,  ze )  =  -^2  f  d2kn(ko,  ,  (45) 

a  2-D  wavenumber  integral  whose  integrand  has  amplitude  k(&o,  k)  and  phase 


y(k,K,k0;rjoe,zo,ze )  =  (k  -  K)-Lo£  +  K(ko,  Ik  -  K\)(zq  +  ze)  . 


4-3.2  Stationary -Phase  Evaluation 

The  shallow-deployment  conditions,  essentially 

e  d-  ZQi/rw  <  1  , 

facilitate  applying  the  stationary-phase  method  to  estimate  L.  First,  the  result  with  K_  =  Kn{6) 
and  L  =  2  (for  Eq.  (35a))  is  found  explicitly.  The  result  with  K_  =  Kn{6  +  7 r)  and  1  =  1  (for  Eq. 
(35b))  is  obtained  by  inspection.  To  simplify  the  process,  we  introduce  the  dimensionless  variables 

/£  =  k/ko  v  =  K/2k0  ip  =  ty/koToe  (46) 

and  the  function  7(2)  =  \/l  —  x2.  With  these,  L  can  be  expressed  (with  the  dependence  on 
everything  but  /i,  u  left  implicit)  as 

L(z)  =  ik0  (^)2/  7(y)cifc0P0^*K)  (47) 

where 


V>(/£5  £)  =  (M  -  2^)*n(^)  4-  7(1^  ”  2z^|)2e  .  (48) 

Under  the  innocuous  assumption  that  many  wavelengths  separate  the  source  from  the  receiver  (i.e., 
that  ^>1)  the  stationary-phase  estimate  is 


L(e)  « iko  (^ — )  J^-yflal)  |det  K(ql,e) 


ko 


V47T  roe 


-1/2 


,ik0  rat_ib{g_,u) 


(49) 


in  which  the  g_  are  the  stationary-phase  points  (i.e.,  points  where  ip  is  real  and  dip / dji  vanishes), 
and  MX(±,v)  =  —  v) /  8  is  the  phase  curvature  tensor.  To  evaluate  L  for  a  particular 


16 


R.F.  Gragg 


v,  one  must  first  find  each  a,  then  evaluate  the  amplitude  factors  7(|ct|)  and  det  u)  and  the 
phase  t/> It  is  a  straightforward,  though  tedious,  calculation  to  do  this  exactly.  We  omit  the 
details  and  note  these  approximations  for  e  <C  1: 

a  =  2 n{6)  +  (1  -  2e2)n{6t)  +  0(e4)  7(|£l)  =  2^1  +  cos  ¥t  +  0(e2) 


|det  MXz.,k)  1/2  =  2e  +  0(e3)  =  k0r0e  +  0(e2) 


2\  ’ 


so  that  the  stationary-phase  estimate  is 

L(u)  = 


kh\ 


ofoe 


7 rr, 


\/l  +  COS  $£  € 


+ikor0£ 


0£ 


No  summation  is  needed  because  there  is  only  one  stationary  point  for  each  u. 

The  above  result  gets  used  as-is  with  l  =  2  in  Eq.  (35a)  and  with  l  =  1  and  6  — >  0  +  tt  in 
Eq.  (35b).  The  factors  needed  are,  respectively, 


kho 


[L(K,TJdti^iz2)]K=K(oj')n(6)  ~  ^f2^  ^  +  C0S e 

rr*,r/  ^0^01 


+«fcoro2 

Vl  -costfi  e-ifcoro1 


(50a) 

(50b) 


'01 


Although  evaluation  with  \K_\  =  K(u>')  is  called  for,  it  is  actually  unnecessary;  the  outcome  depends 
on  ICs  direction  but  not  its  magnitude. 

4-3-3  CSD  Estimate 

With  Eq.  (50)  for  the  L  factors,  with  Eq.  (42)  for  the  remaining  Green’s  functions,  and  with 
Eqs.  (14)  and  (104)  to  simplify  the  angle  and  frequency  integrals,  Eq.  (35)  reduces  to 


m-  .  ieiko(r°2~ro1^  2(fc0hrms)2  z02  .  f,  zi 

W;*0,*1,Z2)  =  -6(f)  47r2ro;— - - -  X  —  Sln  W 


)  W(02) 

in  terms  of  the  azimuth  integral 

W (0e)  =  j>  dOH (9)  x  ^  (y^l  +  cos(0  -  $£)  +  \fl  -  cos (0  -  ^))  . 


(51a) 

(51b) 

(52) 


4.4  Baseband  Summary 


The  full  baseband  result  through  second  order  is  Cr°e°f  +  Cr20f  +  Cr02f,  which  can  be  written  as 
Clef(L01,L02,f;z0,ZliZ2)  =  ^(/)Q(Eoi)-2:0,2i)Q*(ri02,20,2:2)  (53) 


in  which5 


Q(Loe,zo,ze) 


(54) 


5  The  relative  sizes  indicated  below  the  factors  in  Eq.  (54)  have  been  exploited  in  rewriting  the  expression  according 
to  the  pattern  (S1S2  —  ie2.X2.s1  +  iei.XiS2)  =  ($i  +  iei-Xi)(s2  +  ^2X2)*  +  0(e2)  for  se,Xe  ~  1,  ee  <C  1.  Figure  8  shows 
both  H(6)  and  the  geometric  factor  multiplying  it  in  the  integrand  of  Eq.  (52).  That  figure  confirms  that  W(6t)  ~  1. 
It  is  easily  verified  that  87r(/irms/Ao)2  ~  1  for  /o  <  300Hz,  U  <  20m/s. 
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Fig.  1  —  Bragg  scattering  at  surface  point  p  by  a  surface  wave  component  with  wavevector  K_  =  ±Kn(0).  The  parallel 
dashed  lines  representing  the  wave  crests  are  spaced  A  =  2x /K  apart.  The  symmetric  surface  point  s  is  shown. 


4.5  Sidebands,  Fundamental  Integral 

With  the  Green’s  function  of  Eq.  (41),  the  fundamental  sideband  integral  Eq.  (39)  becomes 


I(K,rjof,zo,ze)  =  J  d2k  exp{m(k,K,  ko'^zo^i)} 


*  =  1,2  , 


whose  integrand  has  unit  amplitude  and  phase 

ko',  rjoe,  zq,  zi)  =  k-rjoe  +  K.(k0, 1 k-  \K\)zq  +  k(/cq,  | k+  \K_\)ze 


(55) 


(56) 


for  any  surface  wavevector  K_.  To  estimate  I(K .  •  •  •)  for  a  given  geometry,  we  will  need  to  determine 
the  range  of  K_  in  which  stationary-phase  fc’s  arise  and  what  their  values  are. 

Before  embarking  on  that,  it  seems  advisable  to  review  the  physical  significance  of  the 
stationary-phase  condition.  Stationarity  (i.e.,  ■  ■  -)/dk  =  0)  amounts  to 


Lot  =  * 0 


k-hK 


+  Z£ 


k  +  \K 


K(k0,\k-kK\)  K(ko,\k+  hK\)  ' 


(57) 


Any  k  that  solves  this  for  the  prevailing  K_  can  be  used  to  define  a  point  p  on  the  sea  surface  via 
Loi  =  l IoP  +  Lpi  (see  Fig.  1)  with 


k-\K 


£o E  _ 

zq  K,(ko, \k-  \K_\) 


and 


T-pp_ 


k+  \K 


zi  n(ko,  |^+  \K_\) 


(58) 


In  plane-wave  terms,  these  correspond  to  propagation  to  and  fromp  with  the  3-D  wavevectors 

kto  =(k—  \K)  -  K(ko,  | k  -  \K_\)£z  and  £from  =  (k+  \K)  +  n(ko,  \ k  +  \ K\)ez  .  (59) 

Upon  redirection  at  p,  the  wavevector  changes  by 

A kd=f  fcfrom  -kto  =  K+{K(k0,\k  +  \K\)  +  K(k0,\k-  \K\) }ez  .  (60) 

Because  of  energy  conservation,  we  have 
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|&to|  —  l&froml  —  h  .  (61a) 

Furthermore,  for  K_  =  ±Kn(6),  one  finds  from  Eq.  (60)  that 

Ahn(9  ±  ±tt)  =  0,  \Ahn(0)\  =  K.  (61b) 

Equation  (61)  is  precisely  the  condition  for  the  signal  to  travel  from  the  source  to  receiver  £  by 
first-order  Bragg  scattering  from  the  surface  wave  with  wavevector  K_  at  point  p  on  the  sea  surface. 
If  the  values  of  r^p  and  r_pg  are  swapped,  the  same  surface  wave  scatters  the  signal  from  a  point 
5  located  at  r$s  =  Lp£  and  r ^  =  r$p.  These  points  p,  s  lie  on  the  port  and  starboard  sides  of  the 
“heading”  r^£.  In  Fig.  1,  they  are  the  vertices  of  the  parallelogram  that  lie  on  opposite  sides  of  the 
diagonal  r q#. 

Stationary-phase  estimates  of  I(K_,  •  *  ■)  for  the  wavevectors  K_  =  ±Kn{9)  provide  the  up/down- 
Doppler  sidebands  ±|/|  in  Eq.  (38)  that  stem  from  first-order  Bragg  surface  scattering.  Equation 
(55)  is  evaluated  next,  and  the  result  is  used  to  determine  the  CSD  contribution. 

4.6  Sidebands,  Stationary-Phase  Evaluation 

This  section  develops  the  stationary-phase  estimate  for  the  fundamental  sideband  integral 
J(Jf, Lo£'izQiz£)  'm  Eq.  (55).  The  procedure  is  like  the  one  in  the  baseband  case  but  necessarily 
more  complicated.  Two  separate  depth/range  parameters  are  involved, 

h  =  —  ,  H  =  —  with  Vi  =  max(f£,  Si)  ,  (62) 

roi  1 

and  we  are  still  concerned  with  shallow  deployment,  i/£  <C  1. 

In  terms  of  the  dimensionless  fi,v_  and  ip  from  Eq.  (46),  the  integral  at  issue  has  the  form 

I(v)  =  (^)2/  d2  (63) 

with  the  dependence  on  everything  but  fi  and  y_  left  implicit  again.  For  1,  the  stationary- 

phase  estimate  is 

I(v)  «  ( V  I  det  M(o r,  v) |-i/y  .  (64) 

\A-KrwJ  “ 

As  in  the  baseband  case,  g_  is  a  stationary  point,  ip  the  (scaled)  phase,  and  M_  the  phase  curvature. 
For  the  sideband  case,  this  phase  and  its  stationarity  condition  are 


ipifh  k)  =  n(9i)-jj,  +  Sadfj. 

fJL  —  V 


n($i )  =  Si 


7(l/£-M|)  +  €£ l(\t±  +  ”\)  ’ 


v\)  +  ea(  \n  +  id) 

H  +  v 


(65) 

(66) 


Estimates  are  needed  for  both  sidebands  in  Eq.  (38).  With  v  >  0,  this  means  evaluating  Eq.  (63) 
at  v_  =  +vn{0)  and  at  y_  —  —vn(0).  In  this  report,  we  deal  explicitly  with  the  upper  sideband  (the 
plus  sign),  and  obtain  the  lower-sideband  result  by  inspection.  For  any  given  we  first  need  to 
find  the  stationary  point  by  solving  Eq.  (66)  for  /x  within  the  region  where  ip  is  real- valued.  To  do 
that,  we  adopt  a  pair  of  cartesian  axis  vectors 


(ei,e2)  =  [n{B  -  iyr),n(0))  , 


(67) 
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Fig.  2  —  ^-plane  shown  for  the  upper  sideband  case.  The  two  disks  are  the  regions  \p  —  u\  <  1  and  \p  +  v\  <  1.  Any 
stationary  phase  points  must  lie  in  their  overlap,  where  is  real. 


where  the  first  one  points  along  the  crests  of  the  surface  waves  and  the  second  along  their  propa¬ 
gation  direction  (see  Fig.  1).  The  problem  can  now  be  understood  in  terms  of  the  vectors  y_  =  z/e2, 
7  =  7£i  and  the  two  disks  shown  in  Fig.  2.  These  disks  have  unit  radii,  their  centers  lie  at  the  • 
points  where  p  =  rtz/,  and  their  edges  intersect  at  the  “vertices”  o,  where  pi  =  ±7.  Any  stationary 
points  must  lie  within  their  overlap ,  the  \p±v\  <  1  region  bounded  by  so  that  both  j(\p  +  u\) 

and  j(\p  —  u\)  wdl  remain  real.  The  need  for  such  an  overlap  restricts  v  to  values  below  unity. 
Since  8%  and  e#  are  both  much  smaller  than  unity,  Eq.  (66)  can  only  be  satisfied  when  one  of  the 
denominators  on  its  right-hand  side  is  similarly  small.  This  means  that  the  stationary  points  must 
lie  near  the  boundary  of  the  overlap  region.  In  scaled  variables,  Eq.  (58)  for  scattering  point  p 
becomes 


£op=  H-V.  •,  Ip*  H  +  F- 

z0  7(li£-d)  an  *e  7(l/£  +  d)  ' 


(68) 


If  7(| p  —  u\)  is  small  compared  to  unity,  fi  must  be  close  to  the  lower  edge,  and  rop  is,  therefore, 
large.  On  the  other  hand,  if  +  v\)  Is  small,  //  lies  near  the  upper  edge  and  rp£  is  large.6 

The  simplest  analytical  approach  to  finding  the  stationary  points  is  to  adopt  bipolar  coordinates 
(£,77)  since  they  conform  to  the  shape  of  the  overlap  (see  Appendix  A).  (In  adapting  bipolar 
coordinates  to  the  present  problem,  the  function  j(x)  =  \/l  —  x2  has  to  be  evaluated  so  often  at 
x  =  v  that  7(z/)  is  abbreviated  as  7  wherever  it  occurs.)  The  top  and  bottom  edges  of  the  overlap 

6There  are  solutions  for  which  both  y(\p  —  u\)  and  7(|/x  +  v\)  are  small,  and  p,  therefore,  lies  near  one  of  the  vertices. 
These  are  neglected  here  because  geometric  spreading  with  both  rop  and  rpe  large  makes  a  far  smaller  contribution 
to  the  reverberation  (smaller  by  at  least  30  dB  for  all  the  simulations  to  appear  in  this  report). 
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Fig.  3  —  Upper  edge  of  the  overlap  region — right  half 


are  the  coordinate  lines  77  =  770  and  77  =  27t  —  770,  where 

770  =  arccos(— v)  =  |  +  (  (69) 

(  =  arctan(i//7)  .  (TO) 

(  is  the  acute  angle  in  Fig.  2.  The  angle  770  lies  in  the  second  quadrant  (^7r,7r),  while  2n  -  770 
is  in  the  third  quadrant  (7r,  §7r).  The  interior  of  the  overlap  corresponds  to  -00  <  £  <  +00, 
770  <  7?  <  27r  —  770.  The  edges  are  reached  by  the  limits  77  J.  770  and  77  f  27t  —  770  along  const  ant-£ 
lines.  We  use  w  =  and  write  R(0)  to  indicate  a  counterclockwise  plane  rotation  through  0. 


4.6.I  Stationary  Points 


First,  we  examine  the  case  with  large  ry/js*  and  look  for  stationary  points  near  the  upper 
edge  of  the  overlap.  This  edge  is  the  coordinate  line  77  =  770  on  which  /x  is  a  function  of  w  only: 
\x  =  /xQ(w).  At  each  point  on  it,  the  lower  disk’s  unit  radius  vector  f£0(w)  +  u  is  tangent  to  the  w 
coordinate  line  through  that  point  (because  £,77  are  orthogonal  coordinates). 


We  begin  by  focusing  on  the  right  (£  >  0)  half-plane  (see  Fig.  3).  Starting  at  the  upper  edge, 
we  apply  a  small  positive  shift  67]  to  move  downward  along  the  w  line  to  a  point  ji  —  /xQ  +  <5/x  inside 
the  overlap.  To  first  order,  the  inward  shift  6fZ  is  antiparallel  to  /Xq  +  v.  Throughout  the  right 
half-plane,  Eq.  (A5b)  provides  an  exact  result,  // (7*7,77)  =  [1+  wR(v)][ 1“  w&(rl)]~1Z>  from  which 
the  first-order  inward  shift  is 


x  6  r]  = 


2wdR(r])/dr] 

(I  —  wRi?}))2 


7^77  . 
go 


(71) 


The  derivative  of  the  rotation  is  simply  <9 j?(?7) /<9t7  =  R(tt/ 2  +  77),  and  the  inverse  of  1  —  wR(rj)  is 
given  by  Eq.  (A6).  Since  (cos 770,  sin 770)  =  (—177),  one  has 


/x 0  =  (aei  +  27 2we2)/d  ,  /xQ  +  u  =  {aex  +  ce2)/d  ,  =  (a£i  +  ^2)/^  (72) 


6fi=- 


2jw6t] 


{ae^  +  ^2) 


(73) 
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in  terms  of  the  temporary  variables 

a  =  7(1  —  w2)  >0  c  =  is  +  2w  +  isw 2  >  0 

b  —  —v  +  2(1  —  2 v2)w  —  isw2  d  =  1  +  2 isw  +  w2  >  0  . 


(74) 


Variables  a  through  d  have  absolute  values  less  than  4  throughout  0  <  is,  w  <  1,  and  they  are 
interrelated  through  a2  +  c2  —  d2  and  c  —  b  =  2 isd.  Since  \fj^  +  y_ \  —  1  and  \6fi\  1,  Eq.  (68) 
becomes 


—  - - - (1  +  0(|%|))  and  ^  =  1;(1  +  0(|%|)), 


*0  2^/^F 


Z£  y-2(/£0 +  *')•% 

which,  with  Eqs.  (72)  and  (73),  is  equivalent  to 

rtK^m{m+CS2)  and 

The  stationarity  relation  Eq.  (66)  then  reduces  to 

1 


(ae_i  +  &e2)  . 


n(0e)  = 


2y/j wd  l^/2'yis 


^  (ae-L  +  be2)  +  -^=(aex  +  ce2) 


(75) 


(76) 


(77) 


In  cases  of  practical  importance,  the  solutions  for  w  that  emerge  from  this  equation  are  not  ar¬ 
bitrarily  small.  This  can  be  seen  in  the  following  way.  Begin  with  the  definition  v  —  K/2ko, 
invoke  the  gravity- wave  dispersion  relation  K  =  u2 / g  and  take  u  to  be  the  Pierson-Moskowitz 
dominant  frequency  Eq.  (103).  Then  v  has  the  value  is ^  =  u2/2gko  =  \ffij2 0  x  gc^f^f^U2). 
For  5m/s  <  U  <  20m/s,  50Hz  <  /o  <  200Hz,  and  is  =  is the  factor  \f  \J2^v  in  Eq.  (77)  is  roughly 
4  (±3).  This  precludes  w  ~  is2  because  values  that  small  would  leave  the  Si  term  of  order  unity 
and  the  €£  term  much  greater  than  unity,  making  it  impossible  for  n{$i)  to  have  unit  length.  There 
may  be  solutions  with  w  ~  is£,  but  there  are  certainly  none  with  w 

If  we  neglect  the  terms  in  Eq.  (77)  containing  Si/ y/2jis  relative  to  the  much  larger  terms 
containing  e^/ and  form  the  components  in  the  e2  and  e1  directions,  the  outcome  is 


^Aywd  cos  $£ 


J±_ 

V^rj 


c  and  yf^ywd  sin 


*1 


(78) 


Since  a  and  c  are  non-negative,  a  solution  for  w  and  8rj  exists  only  when  the  angle  $1  lies  in  the 
first  quadrant  (0,  ^7r).  With  6g  eliminated,  Eq.  (78)  reduces  to  acos^  —  csin-tf^  =  0,  a  quadratic 
equation  in  w  whose  physical  solution  is  readily  found  to  be 


w 


7  —  sin  $£ 


7  cos  •di  +  is  sin  $£  5 


(79) 


provided  the  angle  is  confined  to  0  <  ^  <  r  —  £.  The  angular  restriction  ensures  that  cos'd  £  >  v 
and,  ultimately,  that  w  lies  in  the  physical  range  0  <  w  <  1.  That  entire  range  appears  accessible 
only  because  the  unattainably  small  values  w  ~  e2  have  not  yet  been  excluded. 

Squaring  and  adding  the  two  components  of  Eq.  (78)  produces  6g  =  (ej/iy)  X  (d/w)  which, 
with  w  from  Eq.  (79),  becomes 
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Fig.  4  —  Admitted-angle  sector  for  the  upper  edge,  for  Vi  =  0.1.  Curve  (A)  is  $£  =  arccosz/.  The  area  above  is  the 
nonphysical  w  <  0  region.  Curves  (B,C)  are  the  exact  and  approximate  forms  of  de  =  Ht{v)  from  Eqs.  (81)  and  (85). 
The  maximum  wavenumber  v  —  1  —  vt  is  indicated  by  the  •.  The  dashed  line  shows  the  range  of  admitted  $£  values 
at  v  =  0.2  (i.e.,  0  <  dt  <  1.25  w  70°). 


Because  the  azimuth  has  been  restricted  so  that  cos#*  —  v  >  0  (i.e.,  #*  <  arccos(zz)  =  —  C)? 

we  know  that  Srj  is  positive.  Further  restriction  is  needed  to  guarantee  that  it  is  small.  Suppose 
we  want  to  ensure  that  8rj  <  (ej/ve)  x  t/2;  we  will  need  to  require  cos$£  —  v  >  V£,  which  means 
reducing  the  maximum  allowable  $£  from  ^7r  —  (  to 

Hl(v)  =  arccos(z^  +  Vjj)  .  (81) 

This  raises  the  minimum  value  of  w  obtainable  via  Eq.  (79)  from  zero  to  z^/272 — the  w  ~  v# 
minimum  that  was  anticipated — and  lowers  the  maximum  v  from  1  to  1  -V£,  For  a  given  geometry, 
stationary-phase  points  of  the  present  type  occur  only  for  surface  waves  in  the  sector  0£  <  6  < 
0£  +  H£{v).  Their  bipolar  coordinates  are  given  by  Eqs.  (79)  and  (80). 

The  rest  of  the  upper-edge  stationary  points — those  in  the  left  half-plane — are  found  in  the 
same  way,  and  the  lower-edge  points  follow  suit  with  only  minor  differences.  We  omit  the  details 
and  just  present  the  results. 

In  summary,  we  find  the  following  for  the  sideband  stationary  points.  All  of  them  have 

7  -  |  sin#* | 


w  = 


7I  cos  #*|  +  v\ sin#*| 
The  upper-edge  stationary  points  have  77  =  770  +  Srj  with 

_  7^1/2 

|  cos$*|  —  V 


(82) 


8rj 


(83) 
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Fig.  5  —  u-0  polar  plots  of  the  forward  admitted  sectors,  for  v  <  1,  with  =  10°,  62  =  50°,  V\  =  0.10,  and 
V2  =  0.15.  Left:  F\ {v)  and  F^iy).  Fz(v)  is  centered  on  Be  and  vanishes  for  v  >  1  —  U£.  Right:  The  intersection 
F{y)  =  Fi(u)  H  F2 (1/),  which  vanishes  at  v  —  umax  (dashed  arc). 


and  correspond  to  near -source  scattering  from  surface  waves  with  6  in  the  azimuth  sector 

Fe(u)  =f  [0e  -  He(u),  0t  +  He(u)}  (84) 

(left  side  of  Fig.  5),  which  is  centered  on  the  forward  direction  0  =  0^  and  has  a  half- width  of 

He( u)  «  \k  -  ((  +  vtji)  .  (85) 

(Equation  (85)  is  an  order-z^  approximation  to  Eq.  (81).  See  Fig.  4.)  The  coordinate  bounds  are 
z^/272  <  w  and  6r)  <  (ej/ug)  x  7/2,  and  the  wavenumber  range  is  0  <  v  <  1  —  vg.  The  lower-edge 
points  have  rj  =  2tt  —  (770  +  6rj)  with 


I  COS  Vg\  —  V 

and  correspond  to  scattering  near  the  receivers  and  to  the  surface-wave  azimuth  sector 

6  e  Ae( v)  d=  [7 r  +  6e-  +  6e  +  Ht{v)\  , 


(86) 


(87) 


which  has  the  same  half- width  Hg(u)  but  is  centered  on  the  aft  direction  6  =  7r  +  8g.  In  this  case, 
8rj  <  ( 'Pi/vt )  x  7/2,  but  the  w  bound  and  v  range  are  unaltered. 

4-6.2  Phase  Factors 

The  phase  of  the  fundamental  sideband  integral  is  given  by  Eq.  (65).  To  lowest  order,  the  small 
eg  and  6g  terms  may  be  neglected,  so  that  'ip(a^u)  «  1  —  v\  cos(0  —  0g)\  for  surface  waves  in  either 
of  the  admitted  sectors.  Since  cos(0  —  0g)  is  positive/negative  in  the  forward/aft  sectors, 


V’(m^)  = 


1  —  ucos(0  —  0i) 
1  +  2/ cos  (6  —  6e) 


6  6  Fe(u) 

0  e  Ae(u)  . 


(88) 
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4-6.3  Amplitude  Factors 

Straightforward  differentiation  produces  the  general  result 

M(  ^  ft  /1  ,  f -1  I  (g  +  g)(g  +  g)’j 

—  7(|£-d)l=  72(k“d)  i  7(k  +  £|)h  72(I<I  +  r|)  / 


(89) 


For  upper-edge  stationary  points,  we  retain  the  term  containing  l/7(|/i  +  j/|)  1  and  drop  the  one 

containing  1/7(|^  —  v\)  1.  For  lower-edge  points,  we  do  the  opposite.  This  leads  to 


det  M 


{ 


upper  edge 
lower  edge 


(90) 


4.6.4  Synthesis 

With  the  amplitude  and  phase  factors  obtained  above,  the  stationary-phase  value  of  the  fun¬ 
damental  sideband  integral  I(u)  in  Eq.  (64)  is 


l(+2kom(0),!j)6zo,zt) 


*0  Jkor0£  y 

4ttt2 


( 


Zie-iko’/Zoe"n(d) 

z0e+ikol/Z>e'!l.(d) 


■■■  ee  Fe(v) 
■  •  •  0  G  At(v) 


(91) 


4-6.5  Bragg-Only  Constraint 

For  v  <  1  —  Hi(v)  is  positive.  Thus  the  source  signal  can  be  Bragg-scattered  into  the  It h 
receiver  by  surface  waves  with  6  in  either  Fp(v)  or  Ap(u).  Requiring  v  <  1  —  v\  and  v  <  1  —  i>2 
would  allow  both  receivers  to  be  ensonified  in  this  way,  but  that  is  not  enough.  The  surface 
wave  components  doing  the  scattering  are  delta-correlated  in  azimuth  (p.  8),  which  means  that  to 
produce  a  nonvanishing  result,  the  same  component  must  scatter  to  both  receivers.  Since  Hp(v) 
decreases  with  v.  we  will  need  to  reduce  the  maximum  v  to  ensure  this.  We  consider  the  forward 
sector  explicitly. 

For  upper-edge  stationary  points,  we  must  admit  only  surface  waves  that  have  6  in  the  common 

forward  sector  F(v)  d=  F\(y)  fl  FfcO').  Since  Hp{v)  vanishes  at  u  =  1  —  vp,  we  must  require  v  < 
1— max(z/i,  vi)  to  prevent  the  individual  sectors  F\{y),  Fziv)  from  disappearing.  To  ensure  that  they 
overlap,  we  need  |(?i2|  <  H\(y) + H2(v) .  At  o  =  1  — max(i/i,  z/2),  we  have7  H\(v)  +  H2(v)  ~  \/2|i^i2|, 
with  the  following  consequences  (see  Fig.  6).  For  \0\2\  <  \/2l,y12|  (e-g-j  lower  •),  no  separate 
condition  is  needed.  For  |6>i2 1  >  (e.g.,  upper  •),  v  must  be  smaller  than  the  solution  of 

|012|  =  arccos(ri  +  v{)  +  arccos(^  +  02).  That  solution  is8 

VQ  =  y/q  +  Pl2  ~Pl2  ,  (92) 

where  q  =  (1  +  2zqz/2  cos  |0i2|  -  cos2  \0\2\  -  i/2  -  z/f)/2(l  -  cos  |0i2|)  >  0.  Thus,  the  maximum  v  is 


Vmax 


{1  -  max(zq ,  U2) 

vo 


I#12|  < 
1^12!  > 


(93) 


7 1/12  and  U12  are  the  mean  and  difference  of  v\  and  v2. 

8 Though  it  fails  for  \012\  ~  7r,  this  solution  is  valid  for  |#i2|  <  7r/2,  the  range  shown  in  Fig.  6,  which  easily  covers  all 
cases  of  practical  interest. 
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Fig.  6  —  The  Bragg-only  constraint  for  the  upper  edge,  with  v\  =0.1  and  v 2  =  0.15.  Curve  (A)  is  the  larger 
of  H\(v)  and  curve(B)  is  the  smaller,  and  curve  (C)  is  their  sum,  H\ (V)  +  The  □  and  A  indicate 

v  =  l—max(Vi,  V2)  and  v  =  1— min(i/i,  1/2),  respectively.  The  ordinate  of  the  point  o  is  2[^i 2  f  «  0.316,  corresponding 
to  |^i2 1  ^  18°.  The  •  points  are  discussed  in  the  text. 

For  0  <  v  <  vmax ,  we  have  the  nonvanishing  forward-sector  intersection, 

F(y)  =  [-5^  +  (  +  max(6>i  +  v>i/y,02  +  V2/1),  ~  C  +  min(0i  -  17/7,02  -  ^2/7)]  (94) 

(see  the  right  side  of  Fig.  5).  For  the  lower-edge  stationary  points,  the  same  sort  of  analysis  applies 

de  f 

and  the  surface  wave  azimuth  6  is  restricted  to  the  aft  sector  A{v)  =  A\(v)  fl  A2{v)  given  by 

A(v)  =  [\-k  +  (  +  max(0!  +  1/1/7,  #2  +  ^2/7),  |tt  -  <  +  min(0i  -  1/1/7,  #2  -  ^2/7)]  ,  (95) 

which  is  simply  F{y)  rotated  by  7r. 

4.7  Sidebands,  CSD  Estimate 

The  final  upper-sideband  CSD  result  comes  from  Eq.  (38).  It  is 

cS0)(roi,ro2,+|a;|;^,*i,^)  =  S(M)  J dOH(e)eik^^ 

x  r(+2koi/n(e),rQl;z0,z1)  x  I(+2kovn(e),rj02]  z0,  z2)  ,  (96) 

with  the  stationary-phase  estimates  for  the  fundamental  integral  I  taken  from  Eq.  (91).  This  is 
evaluated  over  the  appropriate  range  of  v  for  the  geometry,  with  the  integration  including  the 
admitted  azimuths  for  each  v.  We  have  seen  that  this  reduces  to 

CS°)(E01,Z02,+I"l;*0,*i,*a)  =  Aeik^~r^\B  [  +D  [  1  H(0)d0  ,  (97) 

{  ■'F(v)  J  A(y)  ) 
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(98) 
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This  applies  within  the  range  0  <  v  <  umax,  where  the  source  is  able  to  ensonify  both  receivers  by 
Bragg  scattering  from  the  same  surface-wave  component.  The  admitted  sectors  F(v),  A(u)  contain 
the  azimuths  9  of  all  such  components  traveling  in  the  forward  and  aft  directions,  respectively. 
With  the  relative  receiver  displacement  expressed  as  r12  =  QTlW)?  the  phase  term  in  the  above 
CSD  integral  is  simply  2kour12-n(9)  =  2k0i'Qcos(9  -  <p).  A  change  of  variables  reduces  Eq.  (97)  to 


C^f\rol,rJ02,+Hz0,z1,z2)  =  Aeik^~r^x 

j  ^  d0  { BH{9 )  +  De-i2koS,/cos{d-v)H(6  +  tt)}  ,  (99) 

where  ko  and  K  are  wavenumbers  for  the  acoustic  signal  and  surface  waves  and  v  =  K/2k,Q.  The  sea 
surface  is  described  by  its  dispersion  relation  |o>|  =  0(A),  power  spectrum  $(M),  and  directional 
spectrum  H{9).  The  lower  sideband  at  -\u\  is  exactly  the  same  except  that  H(6)  and  H(9  +  %) 
are  interchanged. 

4.8  Sidebands,  Colocated  Receivers:  the  PSD 

When  the  receivers  are  merged  at  a  common  location  A,  the  CSD  reduces  to  the  PSD, 


Crrf0)(lo*,±|w|;2o,^)  =  C'rtf°)fa>2>toi>±M  ',zo,zi,z2) 


^i=ro2=ro«;2i=*2=z. 


with 

A  =  4>(|a;|)(fco/4vrro*)2,  B  =  (z*/ro*)2,  D  =  (^oAo*)2  • 

The  upper-sideband  PSD  is  then 

c£f\rJO*,+Hzo,z*)  =  A[  d0{BH(9)  +  DH(9  +  tt)}  •  •  •  0  <  v  <  vmax  ,  (100) 

JF{u) 

and  the  lower-sideband  result  Cj2/^(r0*,  —  |w|;  zq,  z*)  is  produced  from  it  via  H (9)  ^  H(6  +  ix). 

4.9  Sidebands,  Symmetries 

In  general,  the  upper  and  lower  CSD  sideband  contributions  are  unequal.  However,  as  Harper 
and  Labianca  originally  noted  for  the  PSD  in  the  bistatic  limit  [3—5],  identical  sidebands  do  result 
in  one  special  geometry — when  the  receivers  lie  in  the  crosswind  direction  from  the  source  (i.e.,  in 
the  x  =  0  plane9).  It  is  not  difficult  to  show  that,  for  receivers  anywhere  in  that  vertical  plane,  the 
interchange  H{9)  ^  H(9  +  7r)  leaves  the  F  and  A  integrals  in  Eq.  (97)  unaltered.10  Thus,  in  that 
one  particular  configuration,11  the  sidebands  are  identical  in  both  amplitude  and  phase. 

A  different  kind  of  symmetry  arises  from  the  transformation 

9l,92  - ►  7T  -  01, 7T  -  92  ,  (101) 

which  reflects  rj  and  r2  across  the  x  (i.e.,  crosswind)  axis.  With  the  CSD  expressions  obtained 
here,  only  trivial  manipulation  is  needed  to  reveal  that  the  effect  of  Eq.  (101)  is  simply  to  swap  the 
sidebands.10 

9 References  3  and  5  cite  the  *  =  0  plane;  Ref.  4  cites  the  y  =  0  plane,  but  this  is  presumably  just  a  typo. 

10This  also  requires  the  fairly  innocuous  assumption  of  crosswind  symmetry,  H (0)  =  H(—8). 

11  Of  course,  this  also  occurs  trivially  in  any  configuration  with  a  completely  isotropic  sea  surface  fj,  =  0. 
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The  PSD  sidebands  are  also  unequal  except  in  certain  special  source/receiver  configurations. 
Obviously,  one  such  configuration  is  the  bistatic  limit  (g  -+  0)  of  any  configuration  with  equal  CSD 
sidebands  (e.g.,  any  receiver  located  in  the  x  =  0  vertical  plane).  The  PSD  also  has  symmetric 
sidebands  whenever  the  source  and  receiver  are  at  the  same  depth  (since  zq  =  z*  implies  B  =  D), 
a  fact  recognized  early  on  by  Harper  and  Labianca  [2-6]. 

5.  SIDEBAND  SIMULATIONS 

This  section  presents  several  CSD  sideband  computations  for  an  ocean  of  uniform  sound  speed, 
Co  =  1500m/s.  The  sea  surface  is  characterized  by  expressions  for  the  power  and  directional  spectra 
<&(u)  and  H(9)  that  have  found  theoretical  or  experimental  support.  The  power  spectrum  is  the 
archetypical  one  for  fully  developed  wind-driven  seas— the  Pierson-Moskowitz  distribution  [26], 

$(<*>)  =  ~rexP (~f (w<*/w)4)  •••  w>0  (102) 

(see  Fig.  7),  which  modifies  Phillips’s  a;-5  behavior  [27]  by  introducing  a  dominant  frequency 

vd  =  (i0)1,4g/u .  (103) 

( U  is  the  wind  speed  at  19.5  m,  g  is  the  gravity  constant,  a  =  0.0081  is  the  Phillips  constant,  and 
/?  =  0.74  is  a  second  empirical  constant.)  The  associated  mean-square  waveheight  is 

=  =  (104) 

The  directional  spectrum  is  an  empirical  “cosine-power”  form  inferred  from  at-sea  data  by  Longuet- 
Higgins  and  coworkers  [28,29]  (see  Fig.  8), 

H{9)  oc  |cos(0/2)|2ai  .  (105) 

The  normalization  is  §d6H{6)  =  1,  and  the  directionality  index  //  has  a  non-negative  frequency- 
dependent  value,  typically  about  4  but  sometimes  less  than  1.  Unlike  such  forms  as  H{6)  oc  cos2(0) 
[12],  this  is  not  upwind/downwind  symmetric  (except  in  the  isotropic  pi  =  0  case). 

The  simulations  explore  the  influence  of  source/receiver  geometry  for  various  {7,  /o,  and  pi.  In 
each  case,  values  are  chosen  for  these  three  parameters  and  results  are  then  generated  simultaneously 
for  an  entire  series  of  receiver  geometries.  (The  computation  is  vectorized  over  geometry.)  For 
each  geometry,  the  upper  and  lower  sidebands  are  evaluated  at  80  frequency  points  each,  and  the 
numerical  azimuth  integration  is  done  by  the  trapezoid  method  with  2°  step  size.  Simulations  are 

identified  by  figure  number  (e.g.,  the  CSD  from  simulation  12  is  in  Fig.  12).  These  figures  all  use 

a  common  format.  The  upper  half  is  a  plot  of  the  top12  30  dB  of  the  magnitude  of  the  CSD,  with 
decibel  values  indicated  in  a  color  scale  on  the  right.  Overall  levels  are  low  because  source-to-surface 
and  surface-to-receiver  propagation  losses  have  not  been  removed,  as  is  done  in  computing  surface 
scattering  strengths.  However,  all  the  computations  are  done  at  the  same  average  source-receiver 
range  (5  km),  so  the  results  can  be  directly  compared  with  each  other.  The  lower  half  of  each 
figure  shows  the  CSD  phase  in  degrees,  using  a  360°-periodic  color  scale.  The  ordinate  of  each 
plot  simply  enumerates  the  geometries  used.  The  abscissa  is  the  Doppler  frequency  shift.  Its  axis 

extends  in  both  directions  to  a  maximum  value  fmax  d=  £}(2ko)/2ir  =  y/ g  fo/ ft  cq,  corresponding  to 
12  Values  more  than  30  dB  below  the  maximum  are  not  rendered.  They  appear  as  patches  of  white  space. 
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Fig.  7  —  The  Pierson-Moskowitz  power  spectral  density  <$(u;),  Eq.  (102),  with  the  parametric  dependence  on  wind 
speed  U  shown.  Contours  are  at  10-dB  intervals.  The  dominant  frequency  Eq.  (103),  is  also  shown  running  along 
the  ridge.  For  this  plot,  the  largest  w  corresponds  to  0.5  Hz  and  the  smallest  U  used  is  3  m/s. 


Fig.  8  —  Solid:  the  Longuet-Higgins  “cosine-power”  directional  spectrum  H (0),  Eq.  (105),  for  integer  values  of  from 
0  (circular)  to  5  (most  elongated).  ^  «  4  is  typical.  The  wind  blows  left-to-right,  so  0°  is  the  downwind  direction. 
The  dashed  line  shows  a  factor  in  Eq.  (52). 
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as  described  in  the  text 


the  v  =  1  cutoff  of  Bragg  scattering.  Although  there  is  no  minimum  Doppler  shift  in  the  theory,  a 
de  facto  fm{n  occurs  where  the  30  dB  window  used  in  rendering  these  figures  cuts  off  the  exponential 
low-frequency  rolloff  of  the  surface  power  spectrum.  These  figures  are  all  blank  in  the  baseband 
\f\  <  fmin •  If  should  be  remembered  that  this  is  not  because  the  zero-Doppler  contributions  are 
absent,  but  because  only  the  sideband  contributions  are  being  plotted. 

The  wind  speeds  and  frequencies  used  are  represented  by  the  five  circles  in  Fig.  9.  With  one 
exception,  these  are  chosen  to  lie  between  the  solid  and  dashed  lines  in  that  figure.  The  dashed 
line  is  taken  from  Ogden  and  Erskine  [1,  Fig.  21]  and  represents  their  estimate  for  the  onset  of 
significant  whitecap  formation.  We  stay  below  it  so  that  surface  reverberation  can  reasonably  be 
attributed  to  interface  scattering.  The  solid  line  represents  v&  —  1.  Since  V&  varies  inversely  with 
/o  and  Uy  one  has  <1  only  above  the  line.  Bragg  scattering  from  a  surface  component  requires 
v  <  1,  so  the  region  below  the  solid  line  should  yield  much  lower  reverberation  levels  because  the 
dominant  surface  waves  cannot  contribute  there. 

The  impact  of  \i  is  foreshadowed  in  Fig.  10.  This  figure  presents  a  plan  view  of  the  surface 
points  (*)  where  dominant  waves  scatter  the  signal  to  a  typical  receiver.  Those  on  the  left  side 
lie  relatively  close  to  the  source  and  seem  to  form  one  branch  of  a  hyperbola.  The  ones  on  the 
right  are  grouped  around  the  receiver  and  appear  to  coincide  with  the  other  branch.  The  line 
segments  indicate  directions  of  travel  for  the  pair  of  counterpropagating  waves  that  participate  in 
the  scattering  at  each  *.  Their  lengths  are  H(6)  and  H(6  +  7r),  where  6  and  6  +  tv  are  the  headings 
of  the  two  surface  waves.  Different  colors  indicate  contributions  to  the  upper  and  lower  sidebands 
(USB,  LSB). 

5.1  Horizontal  Geometries 

In  each  of  the  first  ten  simulations,  a  horizontal  reference  point  =  (rc,0c)  is  chosen  and 
the  receivers  are  laid  out  at  a  common  depth  in  the  following  way.  The  first  placement  (with 
geometry  index  1)  has  r^  and  r2  closely  bracketing  r^.  As  the  geometry  index  increases,  Zi  and  r2 
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Fig.  10  —  Plan  view  of  the  surface  scattering  points  for  the  v  =  Ud  components  of  a  Pierson-Moskowitz  surface 
spectrum.  The  azimuth  is  sampled  in  2°  steps,  /o  =  100  Hz,  U  =  5  m/s,  8(  =  30°,  z0  =  110  m,  Z£  =  90  m,  and 
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separate  uniformly  along  a  straight  line,  keeping  r ^  between  them.  The  result  is  a  series  of  receiver 
geometries  with  a  common  value  of  (p  and  with  steadily  increasing  £,  as  in  Fig.  11.  Parameter  values 
are  always  chosen  so  that  the  geometry  index  is  numerically  equal  to  the  receiver  separation  q  in 
meters.  The  depths  are  not  varied.  For  all  cases,  rc  =  5  km,  6C  =  30°,  g  =  (1  m,  2  m, . . . ,  200  m), 
zi  =  Z2  =  100  m,  and  zq  =  90  m.  Table  1  gives  the  values  of  the  remaining  variables. 


Fig.  11  —  Geometries  for  a  horizontal  simulation.  In  this  example,  geometry  index  =  1, . . . ,  200 


Table  1  —  Input  Parameters  for 
the  Horizontal  Simulations 


Fig. 
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Simulations  12-14  correspond  to  the  upper  o  in  Fig.  9.  The  wind  speed  is  the  maximum  from 
that  figure  and  the  frequency  is  the  highest  value  consistent  with  pure  interface  scattering  at  that 
wind  speed.  The  amplitudes  for  these  simulations  have  certain  main  features  in  common: 

1.  sideband  peaks  at  fd  =  0.09  Hz,  due  to  scattering  by  the  dominant  surface  waves, 

2.  amplitudes  that  drop  off  by  30  dB  before  reaching  fmax  —  0.45  Hz,  and 

3.  a  de  facto  minimum  /mtn  =  0.05  Hz,  which  is  produced  by  the  30  dB  cutoff. 

There  is,  however,  considerable  detail  that  varies  from  one  simulation  to  another. 

In  simulation  12,  the  amplitude  is  largest  where  g  is  one  meter  (at  geometry  index  1).  At  that 
point,  the  two  sidebands  are  symmetric  about  zero  Doppler.  With  increasing  separation,  the  USB 
remains  virtually  unchanged,  but  the  entire  LSB  structure  narrows  considerably  as  g  approaches 
200  m.  This  contraction  is  accompanied  by  the  emergence  of  a  set  of  secondary  peaks — the  ridges 
on  the  downshift  side.  These  appear  to  grow  at  the  expense  of  the  main  peak,  which  falls  by  5  dB 
while  the  adjacent  secondary  peak  rises  to  within  10  dB  of  it.  The  corresponding  phase  is  shown 
in  the  bottom  half  of  the  figure.  The  USB  has  very  little  structure.  Each  g  produces  the  same 
phase  for  all  Doppler  shifts,  and  this  value  simply  loops  through  about  91/ 4  regular  360°  cycles  as  g 
traverses  its  200  m  range.  This  is  easily  understood  in  terms  of  the  horizontal  source/receiver  layout 
for  this  simulation,  which  is  shown  Fig.  11.  The  receivers  lie  at  a  horizontal  angle  \6C  -  <p\  =  45°, 
with  the  radial  to  the  source.  Given  the  signal’s  15  m  wavelength,  direct  line-of-sight  propagation 
would  yield  arccos(9.25  *  15/200)  =  46.07°,  which  agrees  quite  well  considering  that  it  neglects  the 
deviation  due  to  grazing  the  surface.  The  LSB  phase  is  virtually  the  same  if  one  looks  only  under 
the  main  peak.  The  secondary  peaks,  however,  are  an  altogether  different  matter.  At  any  given  g, 
the  phase  under  these  structures  is  not  even  approximately  independent  of  Doppler  shift. 

For  simulation  13,  (p  is  reduced  to  50°.  The  effect  on  the  CSD  amplitude  is  to  lessen  the 
large- p  narrowing  of  the  LSB  and  to  fill  in  the  troughs  between  its  secondary  peaks,  leaving  only  a 
slightly  wavy  shoulder  on  the  downshifted  side.  The  phase  shows  a  corresponding  smoothing  of  the 
secondary-peak  artifacts  in  its  LSB.  The  USB  now  executes  approximately  121/2  full  phase  cycles, 
corresponding  to  arccos(12.5  *  15/200)  =  20.36°  ~  \0C  —  p\  =  20°. 

Simulation  14  uses  ip  =  0°  (i.e.,  one  receiver  directly  downwind  of  the  other).  The  phase 
executes  about  ll1/2  full  cycles,  which  yields  the  now  familiar  angular  agreement:  arccos(11.5  * 
15/200)  =  30.40°  vs  1 6C  —  ip\  =  30°.  Beyond  that,  there  is  no  more  than  a  slight  quantitative  shift 
to  distinguish  this  from  the  preceding  simulation.  The  sidebands’  amplitudes  are  quite  similar,  but 
the  LSB’s  phase  is  still  Doppler  dependent. 

These  simulations  were  all  rerun,  with  9C ,  4>  — >  n  —  9c.tt  —  <f>.  That  transformation  is  equivalent 
to  Eq.  (101)  and,  thus,  should  simply  swap  the  sidebands  in  each  simulation.  This  is  exactly  what 
happened.  The  results  are  not  included  since  they  are  merely  Figs.  12-14  held  up  to  a  mirror. 

Simulations  15-17  repeat  the  previous  three  simulations  with  U  reduced  to  5  m/s.  This  wind 
speed  corresponds  to  the  lower  o  in  Fig.  9,  and  since  that  symbol  lies  above  the  solid  line,  dominant- 
wave  scattering  should  still  be  a  principal  feature.  Relative  to  simulations  12-14: 

1.  the  levels  are  down  by  about  25  dB,  of  which  15-20  dB  can  be  attributed  to  the  reduction  in 
$(u>d)  as  the  wind  speed  drops  from  15  m/s  to  5  m/s  (see  Fig.  7); 

2.  fmax,  being  independent  of  wind  speed,  is  unchanged; 

3.  fd  oc  U~l,  so  the  main  Doppler  peaks  are  upshifted  a  factor  of  3  to  about  0.27  Hz;  and 

4.  fmin  rises  by  roughly  the  same  factor  to  about  0.16  Hz. 
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The  USB  remains  as  featureless  as  ever,  but  there  is  more  detailed  structure  in  the  LSBs  and 
greater  variation  among  the  simulations. 

In  simulation  15,  the  LSB  amplitude  falls  off  faster  with  increasing  g  than  it  did  in  simulation 
12,  where  the  wind  speed  was  higher.  Indeed,  by  the  time  g  has  reached  20  m,  the  main  peak 
is  already  down  by  roughly  5  dB.  It  is  hard  to  be  exact  about  this  figure,  however,  because  the 
secondary  peaks  cut  across  the  main  peak,  blending  with  it  to  form  a  heavily  corrugated  structure. 
The  LSB  phase  reflects  this  structure  but  otherwise  differs  little  from  simulation  12. 

Simulations  16  and  17  continue  the  trend  noted  in  simulations  13  and  14.  As  (p  is  reduced,  there 
is  a  substantial  fading  of  the  irregularities  in  the  LSB’s  amplitude.  For  g  less  than  a  few  meters, 
the  LSB  is  always  identical  to  the  USB  in  both  amplitude  and  phase.  As  g  increases  further,  the 
main  LSB  Doppler  peak  drops  rapidly  by  about  10  dB  and  then  begins  to  taper  off  more  slowly, 
while  the  USB  undergoes  no  material  change.  Ultimately,  as  < p  reaches  0°  in  simulation  17,  the  LSB 
amplitude  resembles  a  version  of  the  USB  with  the  top  5  dB  of  the  Doppler  peak  removed.  The 
phases,  however,  remain  markedly  different.  At  any  given  receiver  separation  (i.e.,  any  geometry 
index),  the  USB  frequency  components  are  all  in  phase  but  the  LSB  components  are  not. 

Simulations  18-20  correspond  to  the  three  •  symbols  in  Fig.  9.  The  frequencies  fd^fmax ,  and 
the  decibel  level  of  the  CSD  amplitude  are  examined  and  compared  to  extrapolations  from  simu¬ 
lation  16  based  on  the  relations  fd  oc  U-1,  fmax  cx  \//o,  A  ex  /q  and  the  U  dependence  of 
from  Fig.  7.  In  all  three  simulations,  /mzn  remains  approximately  0.27  Hz,  and  the  number  of  phase 
cycles  in  the  USB  remains  accurately  proportional  to  /o. 

Simulation  18,  with  Vd  =  0.50,  lies  well  within  the  Vd  <  1  region. 

1.  The  level  is  —5.5  dB  relative  to  simulation  16.  To  a  good  approximation,  this  is  the  sum  of  the 

—  11  dB  change  in  as  U  decreases  from  5  m/s  to  3  m/s  and  the  +6  dB  change  in  A  due 

to  the  doubling  of  /o. 

2.  The  maximum  Doppler  is  now  0.45  Hz  *  y/2  =  0.64  Hz. 

3.  The  main  Doppler  peak  is  shifted  to  0.27  Hz  *  5/3  =  0.45  Hz. 

The  extrapolations  from  simulation  16  are  all  accurate. 

Simulation  19  is  the  borderline  case,  Vd  =  1.00. 

1.  The  dB  level  change  is  —14,  whereas  the  dominant-wave  prediction  is  —11. 

2.  fmax  =  0.45  Hz,  as  in  simulation  16. 

3.  fd  =  0.45  Hz,  but  the  main  Doppler  peak  is  actually  at  0.40  Hz. 

The  Ud  <  w  part  of  the  surface  spectrum  is  unable  to  Bragg  scatter.  This  reduces  the  dB  level  and 
even  skews  the  main  Doppler  peak  to  a  lower  value. 

Simulation  20  has  Vd  =  2.00. 

1.  The  dB  level  change  is  —33.5,  but  the  dominant- wave  prediction  is 

2.  fmax  =  0.45  Hz  *  Vol  =  0.32  Hz. 

3.  There  is  no  actual  Doppler  peak,  since  fmax  <  fd  =  0.45  Hz. 

This  simulation  is  so  far  below  the  dominant-wave  cutoff  I'd  =  1  that 
the  dominant  wave  picture  of  Bragg  scattering  are  invalid. 

Simulation  21  repeats  simulation  15  with  an  isotropic  sea  surface.  The  sidebands  are  perforce 
symmetric.  The  amplitude  still  exhibits  5  dB  corrugations,  but  the  phase  has  only  mild  ripple 
features  with  none  of  the  wild  behavior  previously  seen  in  the  LSB. 


—11  —  6  =  -17. 


all  extrapolations  based  on 
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Fig.  12  —  CSD  simulation  12 


geometry  index  geometry  index 


39 


44 


R.F.  Gragg 


5.2  Vertical  Geometries 

Here  the  receivers  are  located  at  a  common  horizontal  position  (r^  =  r01  =  £02)  but  different 
depths.  Thus  g  vanishes,  <p  is  irrelevant,  and  we  have 

C^f\rc,+\w\;z0,z1,z2)  =  A  [  d9{BH(9)  +  DH(9  +  tt)}  •••0  <  u  <  umax  (106) 

JF{v) 

A  =  <F(H)  (fc0/47rrc)2  ,  B  =  ( zc/rc )2  ,  D  =  {z0/rc)2  ,  (107) 

where  zc  =  yjz\z^  is  the  geometric  mean  of  the  receiver  depths. 

For  simulations  23-25,  z\  and  z\ 2  vary  from  20  m  to  120  m  in  10  m  steps,  with  z\  >  Z2-  The 
resulting  66  geometries  are  numbered  in  order  of  increasing  zc  (see  Fig.  22).  In  all  three  simulations 
/o  =  100  Hz,  U  —  5  m/s,  and  n  =  4.  (g  =  0  m,  and  <p  is  undefined.)  In  simulations  23,  24,  and 
25,  zq  is  20  m,  70  m,  and  120  m,  respectively.  Since  Eq.  (106)  produces  real-valued  CSDs,  the 
amplitudes  are  plotted  but  not  the  corresponding  phases. 
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Fig.  22  —  Geometry  index  values  for  the  receiver  depths  used  in  simulations  23-25.  Geometries  are  enumerated 
so  that  zc  —  yjz\ Z2  increases  monotonically  with  the  geometry  index. 


Since  the  receivers  are  situated  under  a  highly  directional  field  of  surface  waves  almost  straight 
downwind  from  the  source,  it  is  a  good  approximation  to  neglect  H(0  +  w)  in  favor  of  H{6).  Then 
the  upper  and  lower  sidebands  are  given  approximately  by 


C^f\r^,+\oj\;z0,z1,z2)  «  | 

A 

• 

[  d0H{0)  ] 

If(u)  ) 

x  B 

(108a) 

C<£f\Lc,-\v\;z0,zi,Z2)  «  | 

k 

[  dOH(9)\ 

J 

\  x  D  . 

(108b) 

This  explains  why  the  upper  sideband  is  independent  of  zq  but  increases  with  the  geometry  in¬ 
dex  (i.e.,  with  zc)  in  exactly  the  same  way  in  all  three  simulations,  while  the  lower  sideband  is 
independent  of  geometry  index  but  increases  with  zq. 
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Fig.  23  —  CSD  simulation  23 
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Fig.  24  —  CSD  simulation  24 
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Fig.  25  —  CSD  simulation  25 

6.  SUMMARY  AND  CONCLUSIONS 

The  problem  of  multistatic  reverberation  from  a  sea  surface  has  been  posed  with  a  narrowband 
source  and  a  pair  of  receivers  deployed  beneath  the  surface  of  a  fully  developed  wind-driven  sea.  The 
small- waveheight  approximation  was  used  to  express  the  cross-spectral  density  of  the  reverberation 
as  a  sum  of  baseband  and  sideband  integrals  over  the  mean  sea  surface.  Statistical  stationarity  and 
the  customary  angle /frequency  factorization  ansatz  for  surface- wave  moments  were  used  to  reduce 
these  expressions  to  manageable  form — notably  Eqs.  (29)  and  (30)  for  the  sidebands.  The  further 
assumption  of  a  homogeneous  ocean  and  shallow  source  and  receiver  placement  allowed  these  ex¬ 
pressions  to  be  simplified  into  analytic  forms  containing  only  a  single  (azimuth)  integration.  The 
sideband  result,  Eqs.  (98)  and  (99),  was  embodied  in  a  Matlab  script  with  (a)  the  gravity-wave 
dispersion  relation  for  deep  water,  (b)  the  Pierson-Moskowitz  power  spectrum,  and  (c)  an  empir¬ 
ical  directionality  from  Longuet-Higgins.  This  computer  tool  was  used  to  explore  the  parametric 
dependence  of  the  sideband  structure  on  the  depth  and  orientation  of  the  source  and  receivers,  the 
frequency  and  wind  speed,  and  the  directionality  of  the  surface  waves. 

The  aim  of  the  work  has  been  to  provide  a  rapid,  reliable  means  of  calculating  the  reverberation 
CSD,  particularly  the  Doppler  components,  for  essentially  any  pair  of  receivers.  This  capability 
is  essential  in  designing  and  evaluating  beamformers  for  active  operation  against  moving  targets. 
The  output  of  an  array  (whether  horizontal,  vertical,  or  volumetric)  and,  thus,  the  input  for  its 
associated  beamformer  (whether  conventional  or  nonlinear)  is  the  complex  cross-spectral  density 
matrix — the  matrix  of  cross-spectral  densities  for  the  hydrophones  taken  in  pairs. 

One  should  really  not  expect  to  make  many  sweeping  predictions  about  the  operational  per¬ 
formance  of  beamformers  by  simply  pondering  the  structure  of  the  reverberation  CSD.  It  is  too 
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complicated.  The  amplitudes  and  phases  of  its  sidebands  are  nontrivial  functions  of  frequency,  wind 
speed,  and  surface  directionality  in  addition  to  source  and  receiver  position.  Furthermore,  the  pro¬ 
cess  of  creating  an  effective  beamformer  is  an  art  unto  itself,  and  the  CSD  is  only  the  raw  material, 
not  the  finished  product.  Nevertheless,  some  rudimentary  design  ideas  do  seem  to  emerge  from  the 
CSD  analysis  undertaken  here.  For  instance,  it  is  clear  from  the  vertical  geometry  simulations  that, 
if  the  receivers  lie  downwind  (upwind)  under  a  fairly  directional  sea  surface,  then  operating  with 
the  source  at  a  shallow  depth  should  improve  the  LSB  (USB),  though  this  would  not  work  on  the 
USB  (LSB).  However,  the  most  straightforward  approach  might  be  to  look  for  robust  ways  to  re¬ 
move  the  sidebands.  They  are  produced  by  Bragg  scattering  from  a  set  of  points  on  the  sea  surface 
near  each  receiver  and  at  another  set  near  the  source  (the  *s  in  Fig.  10).  One  obvious  approach 
for  eliminating  the  sidebands  is  to  use  directional  sources  and  receivers.  The  simple  expedient  of 
using  a  receiver  beampattern  with  a  vertical  null  of  about  45°  half-width  should  help  significantly. 
The  source  might  be  a  horizontal  linear  array  laid  out  along  the  receiver  direction.  Conical  beams 
could  be  steered  in  the  usual  way,  skipping  the  beam  whose  intersection  with  the  surface  coincides 
with  the  hyperbola  formed  by  the  near-source  surface  scattering  points. 

There  are  several  directions  in  which  extensions  to  this  work  might  be  pursued.  One  of  these 
is  to  allow  the  source  to  have  nonvanishing  bandwidth.  This  has  been  anticipated  and  would 
only  require  an  additional  frequency  convolution.  Another  possibility  is  to  modify  the  surface 
description.  Different  functional  forms  could  certainly  be  used  for  the  directionality  and  power 
spectrum.  It  would  be  quite  easy,  for  example,  to  replace  the  Pierson-Moskowitz  spectrum  with  the 
Toba  spectrum  [30].  A  further  option  is  to  extend  the  treatment  beyond  the  uniform  ocean  model — 
to  ducted  propagation  in  shallow  water,  for  example.  This  would  enhance  the  operational  relevance 
of  the  effort,  though  at  the  cost  of  a  potentially  very  great  increase  in  analytical  complexity. 
Finally,  one  might  improve  the  treatment  of  the  surface  boundary  condition  (e.g.,  replace  the 
small-waveheight  approximation  with  a  small-slope  condition  [30]).  This  would  mean  a  complete 
reformulation  of  much  of  the  analytic  approach,  and  the  level  of  effort  required  is  unclear. 
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Appendix 

BIPOLAR  COORDINATES 


For  any  positive  7,  bipolar  coordinates  (£,r?)  are  defined  as  follows.  The  £  coordinate  lines  are 
parameterized  by 

(x  —  7coth£)2  +  y2  =  (7/sinh^)2  •••  —  00  <  £  <  +00  ,  (Al) 

which  defines  a  circle  with  radius  7/|sinh£|  and  center  ( x,y )  =  (7coth£,  0).  Circles  with  ^  >  0  lie 
in  the  right  half-plane,  and  those  with  £  <  0  are  their  mirror  images  in  the  left  half-plane.  Both 
approach  the  y  axis  as  £  vanishes  and,  in  the  £  — *•  ±00  limits,  they  shrink  down  to  the  “vertices” 
(. x,y )  =  (±7,0).  The  y  coordinate  lines  are  parameterized  by 

x2  +  (y  ~  7  cot  t?)2  =  (7/sin?7)2  0  <  y  <  2tt  ,  (A2) 


which  corresponds  to  a  circle  with  radius  7/ sin 77  and  center  ( x,y )  =  (0,7coth?7).  The  circles  of 
this  second  family  have  their  centers  on  the  y  axis,  and  all  of  them  pass  through  both  vertices. 
Each  circle  is  split  by  the  x  axis  into  two  arcs.  The  arc  in  the  upper  half-plane  corresponds  to  an  y 
value  in  the  interval  [0,  tt],  and  the  arc  that  completes  the  circle  in  the  lower  half-plane  corresponds 
to  that  same  y  plus  7r. 


The  transformations  to  Cartesian  coordinates  (x,  y)  and  polar  coordinates  (r,  0)  are 

7  7 


x  — 


cosh  £  —  cos  y 


sinh£ 


cosh  £  —  cos  y 


sin  y 


(A3) 


2  cosh  £  +  COS  7] 
cosh  £  —  cos  rj 

These  allow  the  plane  vector  r  =  xex  4-  ye 2  =  r^(0)e1  to  be  expressed  directly  in  terms  of  bipolar- 
coordinate  operations:13 


tan  6  =  r-7 

smh£ 


(A4) 


1  +  wR(—y) 

l-wRi-y)1-1 

£<o 

(A5a) 

1  +  wR(+y) 

*"l  —  wR{+y)^~l 

e>0, 

(A5b) 

=  e  Kl.  The  inverse  operators  are  easily  found  to  be 

1  (  c 

—  (1  —  2 w  cos  77  +  w* 

!)  1  (l  -  WM Tyj)  ■ 

(A6) 

The  vertices  are  the  w  —>  0  limits,  r  =  ±7^.  Equations  (A5)  and  (A6)  are  used  directly  in  the 
stationary-phase  analysis  in  the  body  of  the  report. 

13 Ordering  (i.e.,  the  distinction  between  AB_-1  and  B_~l A)  is  irrelevant  here.  Since  these  operators  are  all  2-D 
rotations,  they  all  commute  with  each  other. 
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Fig.  A1  —  Bipolar  coordinates  for  7  =  0.7  Coordinate  lines  are  shown  for:  £  =  1,  2,  3  (circles  in  x  >  0  half-plane); 
£  —  _1}  —2,  -3  (circles  in  r  <  0  half-plane);  7?  =  0.37T,  0.4tt,  0.57T,  0.6tt,  0.77T  (arcs  in  y  >  0  half-plane);  and 
77  =  1.37T,  1.47T,  1.57T,  1.67T,  1.77T  (arcs  in  y  <  0  half-plane).  The  vertices  are  the  two  points  on  the  x  axis  where 
circles  with  different  77  values  intersect. 


